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This  work  describes  investigations  of  several  possible  methods 
which  may  be  used  to  extrapolate  a finite  segment  of  a time  function 
while  minimizing  the  error  introduced  in  the  process.  In  the  frequency 
domain  the  equivalent  problem  is  how  to  improve  the  spectral  resolution 
while  minimizing  the  uncertainty  introduced  by  the  extrapolation  in  the 
spectral  amplitudes  and  phases. 

Several  methods  were  investigated  here  to  obtain  the  direct  decon- 
volution of  the  Fourier  spectrum.  The  limitations  of  the  methods  were 
investigated  in  the  presence  of  random  noise.  These  methods  included: 

(1)  Matrix  inversion  of  the  convolution  matrix.  A rapid  recursive 
method  was  ultimately  shown  to  be  most  effective  and  permit  inversion  of 
very  large  matrices  of  the  convolution  type  (Toeplitz),  The  method  is 
shown  to  be  extremely  sensitive  to  the  presence  of  random  noise. 

(2)  A stripping  method  (called  "nibbling"  here)  in  which  the 
resolution  function  is  cross  correlated  with  the  observed  spectrum  in 
order  to  find  the  point  of  maximum  correlation.  A small  "nibble"  is  then 
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calculated  (a)  has  the  proper  resolution  function  shape,  and  (b)  is 
centered  at  the  point  of  maximum  correlation.  The  nibble  is  subtracted 
from  the  observed  spectrum  and  the  process  repeated.  The  method  worked 
well  in  the  simple  form  tested  here.  Improved  resolution  may  be  obtain- 
able with  a proposed  extension  of  the  method. 

(3)  The  iterative  technique  suggested  by  Papoulis  was  studied  with 
indifferent  success.  Only  functions  of  time  which  decreased  steadily  in 
amplitude  could  be  extrapolated,  while  sinusoidal  functions  did  not  seem 
amenable  to  this  method.  The  presence  of  random  noise  gave  very  poor 
resul ts. 

(4)  Prolate  spheroidal  functional  decompositions  in  both  continuous 
and  digital  form  were  studied.  Great  difficulties  were  involved  in 
calculating  and  storing  the  eigenvalues  and  eigenfunctions.  Since  the 
extrapolation  method  suggested  by  Slepian  requires  the  use  of  very  high 
order  functions,  this  method  was  found  of  limited  use  for  all  but  very 
low  bandwidth  limit  functions  (very  smooth  functions).  A very  interesting 
proposal  is  made,  however,  for  a frequency  interpolation  method  to  improve 
Fourier  spectral  resolution. 
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CHAPTER  I 
INTRODUCTION 

The  art  and  science  of  deconvolution  is  to  eliminate,  at  least 
partially,  the  finite  v/idth  resolution  function  which  is  convoluted  with 
any  practical  discrete  Fourier  Transform.  This  dissertation  studies  some 
methods  of  resolution  enhancement  (or  time  extrapolation)  of  the  Fourier 
spectrum  of  a detected  signal  [1]. 

Chapter  II  describes  the  spectral  estimation  methods  used  classically 
up  to  1981 . 

Chapter  III  presents  the  direct  deconvolution  (using  matrix  notation) 
to  solve  a system  of  linear  equations  to  recover  the  estimated  frequency 
spectrum.  This  chapter  also  develops  the  use  of  a recursive  technique, 
to  elegantly  achieve  the  deconvolution  of  the  Fourier  spectrum. 

Chapter  IV  deals  with  an  original  method,  called  here  the  "nibbling 
technique,"  to  deconvolve  typical  Fourier  spectra  that  present  peaks  (the 
typical  case  for  sinusoidal  signals). 

In  Chapter  V an  attempt  to  use  an  iterative  technique  to  recover  the 
estimated  Fourier  spectrum  is  investigated. 

Methods  used  to  obtain  the  eigenvalues  and  eigenfunctions  involved  in 
the  wave  equation  written  in  terms  of  prolate  spheroidal  coordinates  are 
dealt  with  in  Chapter  VI,  as  well  as  the  deconvolution  and  extrapolation 
of  the  Fourier  spectrum.  The  deconvolution  matrix  obtained  numerically 
in  Chapter  III  is  here  described  analytically  as  a sum  of  prolate  spheroidal 
functions. 
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In  Chapter  VII  some  thoughts  about  the  use  of  the  prolate  spheroidal 
functions  in  the  frequency  domain  are  considered.  Conclusions  are  the 
subject  of  Chapter  VIII.  Appendix  A presents  the  organigrams  used  in  the 
deconvolution  with  matrix  methods  developed  in  Chapter  III;  while  Appendix 
B presents  the  organigram  for  the  iterative  technique  developed  in  Chapter  V. 
In  Appendix  C a note  about  the  prolate  spheroidal  functions  is  introduced. 

In  Appendix  D a study  of  the  digital  form  of  the  prolate  spheroidal  func- 
tions to  deconvolute  and  extrapolate  a Fourier  spectrum  is  presented. 


CHAPTER  II 

PRESENTATION  OF  THE  PROBLEM  AND  PREVIOUS  WORK 
2.1  Introduction 

Until  the  late  1960's,  a variety  of  different  techniques  were  used 
for  power  spectral  estimation  such  as  the  classical  Blackman-Tuckey 
technique,  which  consists  of  the  following  three  steps  given  a finite 
segment  of  a time  function: 

(i)  estimation  of  the  autocorrelation  function  4>y^(T), 

(ii)  use  of  a smoothing  window  (truncation,  Bartlett,  Hamming, 
Hanning,  etc.)  w(r), 

(iii)  Fourier  Transformation  of  w(x)  <f>y„(T).  This  technique  is  also 
known  as  the  Tapered  and  Transformed  (TT)  technique. 

However,  windowing  or  smoothing  has  the  disadvantage  of  introducing 
ripples  (leakage)  into  the  spectrum  [2],  and  even  if  this  leakage  is 
reduced  by  the  use  of  a properly  designed  window,  it  is  at  the  expense  of 
the  frequency  resolution. 

Two  techniques  to  improve  the  resolution  in  power  spectral  estimation 
were  introduced  just  before  1970;  these  are  the  Maximum  Entropy  Method 
(MEM)  by  Burg  in  1967.  It  also  appeared  in  the  literature  of  geosciences 
by  Lacoss  in  1971,  Ulrych  in  1972,  and  in  mathematical  statistics  by  Parzen 
in  1968  and  1969  [3].  This  treatment  introduces  "the  first  principle"  of 
data  analysis:  "The  result  of  any  transformation  on  the  experimental  data 

shall  incorporate  and  be  consistent  with  relevant  data  and  be  maximally 
non-committal  with  respect  to  unavailable  data"  [2:23]. 
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The  second  technique  referred  to  is  the  Maximum  Likelihood  Method  (MLM), 
which  involves  considering  a classical  optimal  filtering  problem  where  the 
filter  is  designed  to  pass  the  power  in  a narrow  band  about  the  signal 
frequency  of  interest  and  to  reduce  to  a minimum  the  effect  of  interfering 
noise  frequencies. 

A difficulty  that  arises  with  the  maximum  entropy  method  and  the 
maximum  likelihood  method  is  that  one  does  not  know  a priori  if  a sample  of 
a discrete  time  series  corresponds  to  an  autogregressive  (AR)  process  (and 
is  best  treated  by  the  MEM),  or  to  a moving  average  (MA)  process  (that  is 
best  treated  by  a conventional  lag-window  approach  (Blackman-Tuckey) ) , or 
to  a mixed  autoregressive-moving  average  process  (ARMA)  where  an  iterative 
least  squares  algorithm  to  guarantee  a minimum  time  delay  feedback  is  best 
used. 


2.2  Problem  Statement 

A statement  of  the  problem  investigated  in  this  work  is  as  follows. 

(1)  Given  an  infinitely  long  data  set  (fl  bandwidth  limited),  as 
depicted  in  Figure  2.1a,  and  a truncation  function  to  represent  real  (finite 
length)  data,  as  depicted  in  Figure  2.1b,  then  the  observed  data  are  as 
depicted  in  Figure  2.1c. 

(2)  Take  the  Fourier  transform  (F.T.)  to  get  a spectrum.  This 
observed  F.T.  is  not  bandwidth  limited  since  the  sine  function  extends  to 

± <». 

f(t)  X rect(t)  >-F(f)  a sinc(f)  [4] 

convolution  operator. 

This  is  the  usual  situation.  A spectrum  with  a finite  width  resolution 
function  (sinc(f))  is  caused  by  the  finite  length  of  data. 
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Figure  2.1a  bandwidth  limited  function. 


/ 

^ rect  function 

Figure  2.1b  Truncation  window. 


Figure  2.1c  Observed  record. 


Figure  2.1  Obtention  of  a data  record. 
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Let  F(f)  be  the  true  (deconvoluted)  spectrum 
Fobs('f)  be  the  observed  spectrum 


Fobs^^^  = F(f)  8 sine  (f). 


For  n data  points,  this  is  a vector  of  n components:  using  a matrix 

notation,  we  can  write 


where  ^ is  (n  x n)  si nc-f unction  matrix  (convolution  matrix), 

£ is  (n  X 1)  deconvoluted  spectrum  (column  vector),  and 
Fobs  ■'■s  ^ observed  spectrum  (column  vector). 

Therefore,  knowing  F^j^^  and  the  convolution  matrix,  we  may  be  able  to  solve 
for  the  deconvoluted  spectrum  F. 


To  find  out  what  this  convolution  matrix  should  be,  let  us  consider 
the  sampling  of  a signal  x(t)  (using  a sampling  function  samp(t)),  and  its 
Fourier  transform. 


2.3  Derivation  of  the  Convolution  Function  to  be  Used  in 
a Finite  Length,  Discretely  Sampled  Time  Function 


x(t) 


X(f) 


W1 


x(t) 


X(f)  8 F{samp(t)} 


x(t)  ’ samp(t) 


f 

samp(t) 
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The  sampled  function  is 


X5.(t) 


= x(t)  • samp(t) 


usually. 


Xs(t) 


= x(t)  • [ ^ comb  ( ^ ) ] 


(2.1) 


so. 


X-(t)  = I x(nAt)  • 5(t  - nAt) 


n=-o 


and  from  Equation  2.1  we  have 


Xg(f)  = X(f)  B [ comb  (At-f)  ] 


(2.2) 


n=-oo 


At^ 


(2.3) 


A bandlimited  signal  with  Nyquist  limit  f^j  is  obtained  by  passing  the 
signal  through  a low  pass  filter  as  shown  in  Figure  2.2a.  In  the  time 
domain,  this  frequency  limitation  can  be  described  by  the  function  h(t) 

FT 

h(t)'< >-H{f)  as  shown  in  Figure  2.2b. 


y(t)  = Xg(t)  fl  h(t) 
and 

Y(f)  = Xg(f)  . H(f)  . 

The  ideal  low  pass  filter  sketched  in  Figure  2.2b  is 

H(f)  = rect  ( ) . (2.4) 

Using  Equation  2.3  and  the  expression  forH(f),  Y(f)  becomes 
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% 


Figure  2.2a 


n 0 

Ideal  Low  pass  filter 


» f 


X (t) 

^ \ 

h(t) 

X (t)  B h(t) 

^ N. 

x,(f) 

H(f) 

• H(f) 

Figure  2.2b  Single  input-output  system. 


Figure  2.2c  Signal  with  highest  frequency  lower  than  the  cut-off 
frequency  of  the  filter. 


Figure  2.2  Low  pass  filter. 
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v(f)  " [ ^ I X(f  - i )]  • rect  ( 4-  ) 


n=-oo 


At 


2f. 


X(f)  = K J x(f  - i ) • rect  ( 4-  ) 


n=-oo 


At 


2f 


N 


In  this  case,  where  the  highest  frequency  is  lower  than  the  cut-off 
frequency  of  the  filter  (as  depicted  in  Figure  2.2c),  we  have 


and  consequently 


V(f)  = ^ X(f) 


but 


y(t)  = x(t) 


y(t)  = Xg(t)  B h(t) 


Replacing  y(t)  by  its  expression  (Equation  2.6)  we  get 


^ x(t)  = Xg(t)  a h(t) 

and  using  Equations  2.2  and  2.4 

00 

^ x(t)  = [ I x(nAt)  5(t  - nAt)]  a [2f,^  • sine  (2fn^t)] 

n=-oo 

00 

x(t)  = At  x(nAt)  • 2f|^  sine  (2fj^(t  - nAt)) 
n=-oo 

00 

= 2f|^  At  I x(nAt)  • sine  [2fn^(t  - nAt] 
n=-oo 

x(t)  = Xg(t)  a sine  ( ^ ) 

x(t)  = I x(nAt)  • sine  ( ^ ) 

n=-oo 


(2.5) 

(2.6) 


(2.7) 


In  the  case  of  finite  sampling  (Figure  2.3),  the  sampling  function  is 
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samp(t)  = rect(Y)  • [ comb(-^)  ] 

and 

Xg(t)  = x(t)  • sainp(t) 

Xg(t)  = x(t)  • rectCy)  [ ^ • comb(^)  ] 

Taking  the  Fourier  Transform,  we  get 

Xg(f)  = T X(f)  B sinc(Tf)  a comb(Atf) 

= T [ X(f)  a corab(At-f)]  a , sinc(Tf) 

Using  Equation  2.3  above,  wejget 

^s^^^  ^ H ^ X(f  - -^)  ] a sinc(Tf) 

n=-oo 

00 

^s^^^  ^ sinc(Tf)  a I X(f  - . (2.8) 

n=-oo 

This  is  the  Fourier  Transform  of  the  (T)  finite  length  segment  of  X(t), 
sampled  at  At  intervals. 

Let  us  see  what  this  function  looks  like  in  a simple  example: 
Assuming  the  signal  at  hand  to  be  sinusoidal  (Figure  2.4),  we  have 

Xg(t)  = x(t)  • • comb(^) 

00 

= I x(nAt)  5(t  - nAt) 

n=-oo 

and 

Xg(f)  = X(f)  a comb(At-f) 

= H J • 

n — 00 
So 


n 


Figure  2.5  Fourier  Transform  of  the  signal  in  Figure  2.4, 
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x(t)  = sin  t 

Xs'f)  = it  J X(f  - if)  • 

n--oo 

Each  region  (i),  in  Figure  2.5,  has  a contribution  on  the  rest  of  the 
ensemble  of  regions.  Therefore,  the  convolution  matrix  form  which  we  desire 
to  study  should  not  be  obtained  from  a single  or  pure  sinc-function  but 
from  a combination  of  sinc-functions  to  take  care  of  the  contamination  due 
to  the  higher  spectral  orders,  as  sketched  in  Figure  2.6.  The  number  of 
spectral  orders  which  must  be  included  will  vary  with  the  specific  problem. 
Three  cases  should  be  considered. 

(i)  The  bandwidth  limit  is  much  less  than  the  Nyquist  limit.  For 
example,  given  a continuous  data  set  of  finite  length,  let  the 
limit  frequency  be  f^^  = 60  Hz.  If  the  sampling  is  done  at  the 
rate  of  10,000  points/sec  or  a sampling  time  interval  of 

10,000  10»000  Hz 

that  corresponds  to  a Nyquist  limit 

= 5,000  Hz 

h « 

In  this  case  the  convolution  function  ought  to  be  a single  sinc- 
function  or  a sinc-function  with  a small  contamination  from  the 
higher  spectral  orders. 

(ii)  The  frequency  limit  is  equal  to  the  Nyquist  limit.  In  the  case 
of  time  series  with  finite  sampling  where  the  limit  frequency 
is  equal  to  the  Nyquist  frequency  (as  depicted  in  Figure  2.7), 


Nyquist  limits 
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the  resolution  function  should  be  a sinc-function  with  a 
substantial  number  (say  30)  of  spectral  orders.  Here  the 
contamination  from  the  higher  spectral  orders  is  maximum. 

(iii)  The  function  observed  is  (x(t))^  not  x(t).  In  the  optics 
domain,  we  know  that  a film  integrates  over  time  and  is 
sensitive  to  the  intensity,  I,  and  not  the  amplitude  ip. 

I - . 

Here  we  have  an  incoherent  light  source  and  the  resolution 
function  should  be  a pure  sinc-function  rectified  to  only 
positive  values. 

In  the  case  of  a laser  holography,  which  is  a source  of  coherent 
light,  the  resolution  function  should  be  a sinc-function,  since  phase  as 
well  as  amplitude  information  may  be  recorded. 

2.4  Some  Resolution  Improvement  Methods  Used  in  the  Past 
Two  techniques  used  for  power  spectral  density  estimation  are  the 
Maximum  Likelihood  Method  by  J.  Capon  in  1969  and  the  Maximum  Entropy 
Method  by  Burg  in  1967  [5]. 

2.4.1  Maximum  Likelihood  Method 

Assume  the  autocorrelation  function,  p,  of  a random  process  is  known  at 
uniform  time  lags  from  -(N-1)  At  to  (N-1)  At  where  At  is  the  interval  of 
time  (lag  interval). 

P^.  for  i = -(N-1),-(N-2),...,-1,0,+1,...,(N-2),(N-1)  . 

For  a frequency  f(Hz)  and  a lag  At(sec)  we  write 


A = Z^rfAt 
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Also  define  a real,  even  taper  function  by  w,  to  use  in  estimating  spectra. 

With  this  notation,  the  conventional  estimate  of  a power  spectrum  is 

(N-1) 


P = 


1 


ni(N-l)  ' 


-i  Xn 


(i  is  such  that  i^  = -1 ) 


which  is  a discrete  form  of  the  relations  of  the  two  mathematicians,  Wiener 
and  Kinchin,  who  for  the  first  time  called  attention  to  the  equivalence  of 
the  autocorrelation  function  and  the  power  spectral  density  [6]. 

.+  00 

power  spectral  density  $ dx 

aa  xa 

J -00 

autocorrelation  function  <p  (t)  = | $(o))  dco 

XX  ^7T  I 


The  window  function  is  given  by 

(N-1) 

" (2N-1)  ^ 

-(N-1) 


Wn  e 


-iXn 


The  method  proceeds  by  forming  the  autocorrelation  matrix  with  the 

P.  ^ R(N,N) 


Let 


E = d,e^'^,e^'^^ 


(T  = transpose) 


Hence  the  conventional  Bartlett  estimator  for  power  spectra  is  given  by 


Pea  = — R E* 

oa*  m2 


and  the  Maximum  Likelihood  Estimator  is  given  by 


1 


ML  T -1  * 

riL  e'  R I E 


(*  = conjugate) 
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To  show  that  the  Maximum  Likelihood  spectral  estimator  represents 
the  optimum  power  output  from  a filter  that  passes  the  signal  in  a narrow 
band  about  the  frequency  A of  interest  and  rejects  all  other  frequency  com- 
ponents, let  the  discrete  time  series  to  be  filtered  be  y(t)  and  the  output 
of  the  filter  of  length  N be  x(t).  Then  the  sampled  output  is  given  by 

N 

x(KAt)  = I a y((K  + 1 - n)  At) 
n=l 

where  the  a^'s  are  the  filter  coefficients.  Let 
y(KAt)  = + n(KAt) 


where  A is  a constant  and  n is  a random  process  with  zero  mean  independent 
of  A and  A.  The  filter  coefficients,  a^,  are  chosen  such  that  Ae^'^^  is 
passed  while  n(KAt)  is  rejected  in  order  for  the  output  x(t)  to  be  an 
unbiased  estimator  of  Ae^”^^ 


n+1 


Dividing  both  sides  of  the  equality  by  Ae''^*^, 

1 = I a^ 
n=l  " 


I a 


n=l 


or 


1 = n p~''2^  ^-■i(N-l)Aw  J 


1 = (E*)'*'  • a 


The  variance  of  x(KAt)  is 


= E{x(KAt)  - E[x(KAt)]}2 
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N 

= E { J a.  n[(K+1-n)At]}^ 
n=l 

= (a*)"^  R a 


where  R is  the  correlation  matrix  of  the  noise  process.  The  term  is 
minimized  subject  to  the  constraint 

(E*)^  • a = 1 
if 


With  this  expression  for  the  vector  of  the  filter  coefficients,  a, 
the  power  spectral  estimator,  corresponds  to  the  minimum  variance  if 
R is  the  autocorrelation  matrix  of  the  random  process  of  interest  and  the 
frequency  for  which  the  filter  is  designed  is  that  at  which  the  power 
spectral  density  is  estimated  [5], 

2.4.2  Maximum  Entropy  Method 

The  Maximum  Entropy  Method  estimates  the  power  spectral  density  using 
the  N first  lags  of  the  autocorrelation  function. 

Let  r be  the  vector  of  the  coefficients  of  a filter  to  whiten  the 
random  process  under  consideration, 

and  define  P = (p,0,0 Or  such  that 

R r*  = p 

R being  the  autocorrelation  matrix  of  the  noise  process.  The  maximum 
entropy  spectral  estimator  is  given  by 

p = PAt 

ME  ^T  r*  J.T  E* 
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For  a sampled  continuous  random  process  one  assumes  that  all  the 
power  is  in  the  Nyquist  Band  where  = l/2At,  At  is  the  time 

lag  and  1/At  is  the  sampling  frequency.  Again,  the  Wiener-Kinchin 
relation  yielding  the  autocorrelation  function  is 


P = 


PM£(f)  e 


icot 


dt 


and  should  satisfy  the  N constraints 

+f. 


Pn  = 


N 


df  n = 0,1,2,. ..,(N-1). 


-f. 


(X  = 27rfAt) 


The  spectral  estimator  P|^|£(X)  is  chosen  such  that  the  loss  of  entropy 


AH  = 


+f. 


Ln  Pm£(X)  df 


-f, 


N 


IS  maximum. 

If  s(X)  is  the  frequency  response  of  a linear  filter  and 


Pme(^)  = |s(X)|^  , 

AH  will  be  the  loss  of  entropy  resulting  from  passing  a Gaussian  process 
through  the  filter. 

The  choice  of  the  estimator  to  be  maximum  is  equivalent  to  the 
choice  of  a filter  that  could  generate  a process  from  white  noise  which 
would  have  the  observed  autocorrelation  function  and  which  would  have 
maximum  loss  of  entropy,  AH  [7], 


19 


2.4.3  Performance  of  the  Above-mentioned  Power  Spectral  Estimators 
Numerical  experiments  have  been  conducted  by  R.  T.  Lacoss  [7]  using 

exactly  known  autocorrelation  functions.  In  particular,  spectra  for  the 
sum  of  2 independent  narrow-band  processes  have  been  considered.  Each  of 
the  independent  processes  had  a spectrum  given  by 

K 

+ (i2TTf  + a)2|2 

with  a total  power  in  the  lower  frequency  bandpass  process  of  one  unit 
and  that  in  the  other  process  of  2 units;  also  the  peak  frequencies  were 
kept  constant  at  0.15  Hz  and  0.30  Hz  by  a choice  of  the  parameter  6.  He 
found  that  the  peaks  for  the  Bartlett  and  Maximum  Likelihood  spectra  give 
the  correct  relative  power;  but  in  the  case  of  the  Maximum  Entropy  Method 
it  is  the  area  under  the  peak  that  gives  an  indication  of  the  power  as 
shown  in  Figure  2.8. 

2.4.4  Other  Methods 

Taylor  series  expansion.  The  Taylor  Series  expansion  method  proposed 
by  A.  Papoulis  [8]  in  1975  to  approximate  the  Fourier  transform  of  a given 
realization  is  described. 

If  a given  function  f(t)  is  bandlimited,  this  means  that  the  function 
is  entire  and  can  be  expanded  in  a Taylor  series.  Thus 

fit)  = I ^ f^"^  (o)  . (2.9) 

n=-oo  "• 

which  converges  for  any  given  t. 

Assume  that  f^(t)  is  known  in  a finite  segment  about  the 

origin,  and  all  its  derivatives  f^*^^(o)  can  be  found.  This  is  in  prin- 
ciple a method  to  extrapolate  f^(t).  Extrapolation  of  f(t)  is  exactly 
equivalent  to  resolution  improvement  in  F(f).  However,  first  of  all  there 
is  a difficulty  in  evaluating  the  derivatives  about  the  origin  in  the 
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Figure  2.8  Deconvoluted  power  spectra  with  different  windows. 
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presence  of  noise,  and  secondly  a truncation  error  is  introduced  in  the 
computation  of  f(t)  from  Equation  2.9.  To  determine  the  spectrum,  F(w), 
of  f(t)  it  is  not  acceptable  to  tolerate  errors  in  f(t). 

Iterative  technique.  The  iterative  technique  consists  of  these  steps 

(i)  Determine  the  Fourier  transform  F^(w)  of  f^(t). 

(ii)  First  iteration.  Since  the  function  is  known  to  be  bandlimited 
truncate  F(w): 


F-|(w)  = F^(w)  • Pjj(w)  (as  in  Figure  2.9). 

(iii)  Then  find  f-|(t)  = F ^ {F-|(w)}  and  form 

r f^(t),  lt|  ^ T 

= f^(t)  + [f(t)  - f^(t)]  p^(t)  = 

as  in  Figure  2.10,  |t|  > t 

then  find 

Gl(w)=F  {g^(t)> 

(iv),.  At  the  nth  iteration 

then  f^(t)  = F"''  [F^(w)] 

and  g^(t)  = f^(t)  + [f(t)  - f^(t)]  p^(t) 

and  G^(w)  = F {g^Ct)} 


Finally  lim  F (w)  = F(w) 
n^  n 


Prolate  spheroidal  wave  function  expansion.  Another  technique  that 
has  been  introduced  in  the  extrapolation  of  bandlimited  functions  uses 
the  expansion  in  prolate  spheroidal  wave  functions. 
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Figure  2.10  First  step  of  the  iterative  method. 
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These  prolate  spheroidal  functions  are  the  solutions  to  the  integral 
equation  (Homogeneous  Fredholm  Integral  Equation) 


A4>(t)  = 


+ T 


,(.)  ;,x)  dx 


(2.10) 


with  the  following  properties: 

(i)  Eigenvalues.  They  are  real  and  positive.  This  is  expected 
because  the  kernel  of  the  integral  equation  is  symmetric. 


e {0,1} 


1 > X-  > X,  > . . . > X ► 0 

u I n 


n-Kc 


(ii)  Eigenfunctions.  To  a given  eigenvalue  X^  there  corresponds 
an  eigenfunction  4>p(t).  These  eigenfunctions  are  orthogonal 

0,  n m 

(2.11) 


+00 


4>p(t)  (})jj^(t)  dt 


= I 


^ — o 


1 , n = m 


and  (J)^(-t)  = (-1)"  cj)^(t) 


(2.12) 


(iii)  The  eigenfunctions  constitute  a complete  set.  Let  f(t)  be  a 
a-BL  (bandlimited)  function.  It  can  be  expanded  as  follows 


The  coefficients  a^  are  given  by 

+ 00 


f(t)  <i>^(t)  dt 


(2.13) 


(2.14) 


— 00 
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The  integral  Equation  2.10  can  be  rewritten  as  follows: 

A<))(t)  = (j)^(t)  B (2.15) 

f 

convolution  operator 

where  (j>^(t)  = 4>(t)-p^(t)  with  p(t)  being  a truncation  window.  Taking  the 
Fourier  transform  of  Equation  2.15  yields 


X$(w)  = $^(w)-p^(w) 


(2.16) 


Equation  2.16  means  that  the  eigenfunctions  are  bandlimited,  i.e.: 


4>(w)  = <j)(w)-p^(w) 


(2.17) 


Taking  the  inverse  Fourier  transform  of  Equation  2.17  yields 


♦(t)  = .Kt)  a 


or 


4>(t)  = 


+ 00 


(2.18) 


(2.19) 


To  show  the  double  orthogonality  of  the  prolate  spheroidal  functions. 
Equation  2.10  is  written  for  a given  eigenfunction,  <J)|^(t). 

" ♦„(x)  dx  (2.20) 


X ())  (t)  = 
n^n'  ' 


-T 


Multiplying  by  the  operator 

f+oo 

4>|^(t)  { } dt 


yields 
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f + ao 

. + 00 

<J>p(t)  <j)|^(t)  dt  = 

<J)|^(t)  dt 

J —00  . 

-00 

+ T 


sin  g(t  - x) 


J.  /„\  iin  uvt  - X 

- x) 


dx 


-T 


X 6 , 
n nk 


II 

<5  , 

nk 

+ T 


nk 


0,  n k 

1 , n = k 


4>p(x)  dx 


-T 


+ 00 


The  expression  with  the  second  integration  sign  is  just  <j)|^(x)  as  indicated 
by  Equation  2.19.  Hence 


r"*"  T 


0,  n k 


<)>n(x)  4>|^(x)  dx  = ■ 


••-T 


Consequently  there  is  a double  orthogonality  of  the  prolate  spheroidal 
functions. 


In  {-oo,+oo} 

. 

— CO 

(2.21) 

In 

. 

-T 

(2.22) 

A time  limited  function  f^(t)  can  be  expanded  in  terms  of  the  prolate 
spheroidal  functions, 


= I b 4 (t)  p (t) 
^ n=0  ^ 


for  all  t 


with 


+ T 


(2.23) 


(2.24) 


thus  the  function  is  known  at  any  t from  an  .expansion  in  the  interval  -t,t 
alone. 
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On  the  other  hand,  a band! imi ted  function  f(t)  can  be  expanded  as 
given  by  Equation  2.13  above,  and  since  the  eigenfunctions  4>p(t)  are 
orthogonal  and  the  energy  of  the  signal  f(t)  is  limited. 


E = 


f^(t)  dt 


using  the  expression  of  f(t),  the  energy  E is  given  by 
00 


E = 


I ^ *71  I \ 


- 1 a a 

r,  m n in 

n,m 


m 
/»+  00 


4>n(t)  (j)^(t)  dt 


E = 


f^(t)  dt 


II 

6 


nm 


J .00 


E = I af  < « 


n=0 


(finite) 


(2.25) 


From  Equation  2.22,  the  energy  of  f^(t)  (t-timelimited  function)  is 
given  by 

f+  T 

2 00 

f (t)  dt  = I K K ^ - 

~ __n  n n 

J-T  "”0 

00  00 

n=0  ^n  ■finite  then  f(t)  = b^  converges  for  every  t, 

and  the  function  f(t)  so  formed  is  bandlimited. 

00 

However,  since  b^  may  diverge,  it  is  not  always  possible  to  find 
a a-bandlimited  function  f(t)  that  agrees  with  a x-timelimited  function 

yt). 

N 


But 
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is  a bandlimited  function  that  exists  for  all  t in  and  is  a good 

approximation  of  f^(t)  when  N is  large.  Indeed  the  mean-square  error  is 
given  by 


1 im 


•+  T 

[fM(t)  - f (t)]^  dt  = I A b^  = 0 
^ n=N+l  " 


If  the  timelimited  function  f^(t)  is  not  a segment  of  a bandlimited 
function,  then  the  function  fj^(t)  will  diverge  for  |t|  > x since 


lim 

N->oo 


.+  00 
• -00 


f^(t)  dt  = 


00 


The  great  numerical  difficulties  involved  in  generating  these  eigen- 
values and  eigenfunctions  have  restricted  their  investigation.  We  will 
attempt  application  of  the  continuous  and  discrete  prolate  spheroidal 
methods  in  Chapters  Viand  VII  [8,9,10,11 ,12;i3, 14]. 


CHAPTER  III 

DIRECT  DECONVOLUTION  BY  THE  INVERSE  MATRIX  METHOD 
3.1  Development 

One  method  to  deconvolute  (restore,  unravel,  enhance,  de-blur)  a 
signal  from  a record  of  data  representing  a linear  superposition  of  the 
signal  of  interest  and  noise  is  to  use  the  matrix  formulation. 

The  sought  after  spectrum  is  designated  by  "the  true  spectrum,"  or 
"the  estimated  spectrum,"  or  "the  deconvoluted  spectrum"  and  represented 
by  Mg. 

The  mixture  of  signal  and  noise  is  designated  by  "the  observed  spec- 
trum" and  is  represented  by  M^.  The  convolution  matrix  operator  is 
represented  by  A,  such  that 


A • 


The  elements  of  matrix  A are  obtained  from  a sine  function  which  is  desig- 
nated by 


B = sine  (-^) 


J “ ...-2, -1,0, 1,2,, 


The  parameter  c varies  with  the  number  of  points  between  2 consecutive 
nodes.  With  this  formalism  the  deconvoluted  spectrum  is  given  by 


M 


e 


The  procedure  outlined  below  will  be  as  follows: 
(1)  Choose  the  function  to  be  deconvoluted. 


(3.1) 


29 


30 


(2)  Choose  a series  of  standard  Fourier  spectra  to  deconvolve. 

(3)  Examine  the  results. 

This  will  be  done  for  a simple  sine  convolution  function  (corresponding 
to  a finite  length  Fourier  transformed  data  set  with  a very  small  sampling 
interval);  (3.2)  a spectrum  free  of  noise;  and  (3.3)  a noise  contaminated 
spectrum.  A signal  with  two  frequencies  and  a complex  signal  are  the 
subject  of  parts  (3.4)  and  (3.5).  In  (3.6)  we  try  to  reduce  the  noise  in 
the  estimated  spectrum  by  increasing  the  length  of  the  realization.  In 
(3.8)  the  case  of  a sum  of  sine  functions  corresponding  to  a data  set 
which  is  more  coarsely  sampled  is  treated.  In  part  (3.9)  the  convolution 
matrix  is  inverted  using  a Toeplitz  recursive  algorithm. 

3.2  Deconvolution  of  a Spectrum  with  No  Noise  Added 

3.2.1  Example  1 

The  spectrum  to  be  deconvoluted  was  obtained  using  program  SIC2.FTN 
to  generate  a sine  function  with  a maximum  at  point  25  with  a data  record 
of  50  points.  And  using  program  DCD.FTN,  based  on  a Gauss-Jordan  algorithm 
with  complete  pivoting,  the  deconvolution  was  performed  to  obtain  a delta 
function,  6,  at  point  25. 

3.2.2  Example  2 

The  spectrum  to  be  deconvoluted  was  calculated  as  the  superposition 
of  2 sine  functions  with  maxima  located  at  2 distinct  points.  This  was 
performed  using  program  SZ.FTN  to  generate  the  sine  functions  at  points 
15  and  35.  The  result  of  the  operation  A ^ * Mq  realized  using  program 
DCD.FTN  was  as  expected,  2 5 functions  at  points  15  and  35.  The  same 
thing  was  repeated  with  sine  functions  located  at  points  10  and  40  and 
then  points  20  and  30.  The  results  were  again  as  expected,  delta  functions 
at  points  10  and  40,  and  20  and  30,  respectively. 
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3.2.3  Example  3 

At  this  point,  instead  of  considering  a sine  function  with  a maximum 
occurring  at  an  exact  division  such  as  sinc(iTJ/c)  a displacement,  E,  was 
introduced  say  sinc[ir(J  + E)/c]  with'E  = 0.1,  E = 0.2,  E = 0.3,  E = 0.4. 

The  deconvolution  of  such  spectra  gave  delta  functions  with  added  ripples 
over  the  whole  spectrum. 

3.3  Deconvolution  of  a Spectrum  with  Noise  Added 

The  effect  of  noise  on  the  deconvolution  of  a known  spectrum  was 
investigated  here.  At  first,  the  observed  spectrum  was  obtained  as 
follows.  A 50  points  record  was  generated  by  a sine  v/ave  sin(7rJ/5)  to 
which  a 100  zeroes  record  was  added  on  each  side  and  then  the  whole  set 
was  Fourier  Transformed  (F.T.)  to  obtain  a sine  function  with  maximum 
located  at  point  26.  This  will  produce  a spectrum  of  many  orders  and  is 
different  from  previous  parts.  Using  program  DCD.FTN  to  deconvolute  the 
spectrum,  the  result  was  a delta  function  at  point  26.  Gaussian  random 
noise  with  adjustable  standard  deviation  and  mean  was  generated  of  50  points 
length  that  was  added  to  the  sine  wave  with  100  zeroes  added  on  each  end 
of  the  summation  record  and  the  whole  was  F.T.  to  obtain  the  observed 
spectrum. 

Different  cases  were  investigated  among  which 


Case 

1 

mean  = 0 

SD  = 0.10 

Case 

2 

mean  = 0 

SO  = 0.25 

Case 

3 

mean  = 0 

SD  = 0.50 

The  deconvoluted  spectra  were  noise  contaminated.  A delta  function 
located  at  point  26  was  distinct  in  Cases  1 and  2 but  got  buried  in  the 
noise  for  Case  3. 
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3.4  Deconvolution  of  a Signal  with  Two  Frequencies 
Plus  Added  Noise 

An  input  signal,  C,  is  considered  here,  where  C = sin(TrJ/2.5)  + 
sin(TTj/5)  + random  noise  with  100  zeroes  added  at  both  ends.  The  signal, 
C,  was  then  Fourier  Transformed  and  truncated  to  50  points  with  two  peaks 
at  points  15  and  40.  The  spectrum  was  then  deconvoluted  to  obtain  two 
peaks,  distinct,  but  with  about  50%  noise  contaminating  the  spectrum. 

This  was  performed  using  programs  PS.FTN  and  FOU.FTN. 

The  same  thing  was  repeated  with  addition  of  a random  noise  having 
10%  standard  deviation  and  zero  mean  to  signals  A and  B.  As  expected, 
the  peaks  were  in  the  spectrum;  however,  the  noise  was  higher  than  the 
previous  case. 


3.5  Tentative  Conclusions 

Simple,  short  length  Discrete  Fourier  Transformed  (DFT)  spectra  may 
be  deconvoluted  using  matrix  inversion  methods. 

(1)  The  sine  resolution  function  is  converted  to  a delta  function 
if  the  peak  of  the  DFT  is  located  at  one  of  the  DFT  harmonic 
frequency  points. 

(2)  If  the  sine  peak  is  centered  between  the  frequency  points,  im- 
perfect deconvolution  and  added  ripples  are  caused. 

(3)  The  presence  of  random  noise  rapidly  degrades  the  deconvolution 
process. 

(4)  The  inverse  matrix  to  the  sinc(f2t)  convolution  matrix  looks 
encouragingly  like  a convolution  matrix  itself. 

There  may  be  problems  with  computer  storage  capacity  when  handling 
such  matrices.  Fortunately,  the  ^matrix  is  of  persymmetric  Toeplitz 
type  with  the  form 
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^n-1 ®0 

Any  convolution  operation  involves  such  a matrix  and  compressed 
storage  is  possible  (e.g.,  symmetric  mode  matrix  storage).  For  large 
matrices,  the  storage  requirements  may  still  be  formidable.  Section 
(3.9)  shows  how  this  may  be  avoided  using  a recursive  inversion  method. 

It  was  our  intent  to  examine  the  matrix  A“^  to  see  if  it  had  the 
persymmetric  Toeplitz  form.  If  so,  then  it  would  represent  a convolution 
operation.  The  function  g(t)  that  forms  the  basis  of  the  matrix  A"^ 
would  then  be  a deconvolution  function. 

Thus,  while  it  is  well  known  that 


f(x) 


r 

f(t)  sine  (t  - t)  dt 

CO 


is  always  at  least  as  wide  as  the  main  lobe  of  sinc(flt),  nevertheless  a 


function  g(t)  can  be  found  such  that 

• T 

g(t)  sine  n(x  - t)  dt  = 6(x) 

This  "finite  deconvolution  function"  may  then  be  used  to  improve  the 
resolution  of  Fourier  Transformed  spectra. 
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Larger  matrices  must  be  examined  to  correspond  to  real-life  situations 
and  to  reduce  numerical  errors  caused  by  describing  the  functions  involved 
with  a coarse  mesh. 

3.6  Large  Matrix  Inversion 

The  length  of  the  observed  realization  was  increased  from  50  to  250 
data  points. 

To  begin,  a sine  function  was  generated  and  stored  in  vector  B,  250 
points  long,  to  form  the  convolution  matrix,  A(250,250).  To  invert  a 
matrix  of  this  size  was  beyond  the  capacity  of  the  POP  11/34  computer 
previously  used.  At  this  point,  it  was  decided  to  use  the  mainframe 
university  computer  and  store  the  results  in  the  POP  11/34.  Using  an  IMSL* 
subroutine,  LINV2F  based  on  a Gauss-Jordan  algorithm  with  complete  pivoting 
in  double  precision,  the  inverted  convolution  matrix,  A^^,  was  obtained 
and  stored  in  a compressed  vector  form,  in  the  POP  11/34  using  subroutine 
LOC.FTN,  in  (1/2)  N(N+1 ) memory  locations;  or,  since  N is  equal  to  250,  in 
31,375  memory  locations  of  four  bytes  each. 

The  input  observed  spectrum  was  obtained  from  the  Fast  Fourier  Trans- 
form (FFT)  of  a 50  points  long  sine  function  with  100  zeroes  padded  on 
each  side  of  the  record: 

7T  1 

{ 100  zeroes  + sin(-y)  + 100  zeroes  } 

To  perform  the  deconvolution,  the  inverted  matrix,  must  be  read  from 
the  POP  11/34  memory.  But  this  was  not  possible  because  the  maximum  num- 
ber of  data  points  that  could  be  read  into  RAM  memory  was  about  16,000 
while  the  total  number  of  entries  was  31,375. 


*International  Mathematics  Subroutines  Library 
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To  avoid  this  problem  the  following  procedure  was  devised  to  recover 
the  estimated  spectrum.  The  inverted  matrix.  All,  was  read  one  row  at  a 
time  (i.e.,  1st  row,  2nd  row,  etc,)  and  then  the  deconvolution  was  realized 
at  each  step. 

^(250,250)  . Mq(250)  = Mg(250) 

t ^ t 

observed  spectrum  estimated  spectrum 


1=1 


‘12  ‘'13 


®i!250^ 


0,1 

M , 
e,l 

0,2 

^0,249 

^0,250 

Using  the  REWIND  operation,  the  250  points  for  each  row,  I,  were  read 
from  the  top  of  the  stored  vector,  A"\  and  this  was  done  for  the  values 
of  J = 1,2,. ..,250. 

This  was  performed  with  programs  RAF.FTN,  SKIP.FTN,  and  OC.FTN 
according  to  the  organigram  of  Figure  A-1  of  Appendix  A.  The  subroutine, 
LOG,  used  to  determine  a vector  subscript  that  corresponds  to  an  element 
A(I,J)  of  the  symmetric  matrix  A stored  in  a compressed  form  gave  a real 
output,  M. 

However,  to  obtain  the  first  point  of  the  deconvoluted  spectrum  it 
took  about  34  minutes  and  20  seconds.  Assuming  a linear  proportionality. 
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for  the  250  cases  it  would  have  taken  approximately  six  days  to  achieve 
just  one  run  (i.e.,  {I  = 1,2,. ..,250;  J = 1 ,2, . . . ,250}) . 

To  reduce  the  computer  run  time  and  the  rewinding  process,  another 
algorithm  must  be  used.  The  matrix  being  symmetric,  the  determination  of 
the  location  of  the  diagonal  element, for  row  I will  determine  the  other 
elements  of  this  row  in  only  one  sweep  of  the  vector. 
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KZO)  = 1 J = 1,2, 3, 4, 5, 6,.. 

KZ(J)  = 1,3,6,10,15,21 

•Location  of  the  diagonal  elements 

J = 2,250 

KZ(J)  = KZ(J  - 1)  + J 
•Location  of  the  elements  of  row  I 
J = 1,1 

LOC(J)  = KZ(I)  - I + J 
J = (I  + 1),250 
LOC(J)  = L0C(J  - 1)  + (J  - 1) 

•Reading  elements  of  row  I 

J = 1,(L0C(1)  - 1) 
read_  'JUNK' 

J =1,1 
read  B(J) 

J = (I  + 1),250 

r LL  = (L0C(J  - 1)  + 1),(L0C(J)  - 1) 

^ read  'JUNK' 
read  B(J) 


The  organigram  is  in  Figure  A-2  of  Appendix  A. 
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Program  RVl.FTN  was  used  to  implement  the  preceding  scheme,  reading 
one  row  of  the  inverted  matrix  at  a time  and  obtaining  the  correspondinq 
recovered  point.  Running  the  program  for  the  first  row  (i.e,,  I = 1 and 

J 250)  took  about  2 minutes  47  seconds,  and  the  result  is  the 

same  as  the  one  obtained  previously. 

Present  Scheme  Previous  Scheme  Same  Result 

2 min.  47  sec.  34  min.  20  sec.  Mg(l)  = 0.048785 

Now,  running  the  program  for  the  250  points  of  the  recovered  spectrum 
took  approximately  11  hours  and  30  minutes.  Again,  this  was  a very  long 
time  and  another  algorithm  had  to  be  devised. 

The  convolution  matrix.  A,  is  a persymmetric  matrix  and  Toeplitz. 
The  inverse  matrix,  A_|_,  is  also  a persymmetric  matrix  of  the  following 
form: 
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Theory  predicts  that  the  inverse  of  a symmetric  matrix  is  also 
symmetric.  This  had  been  checked  with  a Toeplitz  persymmetric  matrix 
(50  X 50)  and  the  result  was  also  a persymmetric  matrix.  Therefore, 
the  knowledge  of  N(N  + 2)/4  elements  of  the  matrix  determines  all  the 
matrix.  In  the  above  case  A(8,8),  the  underlined  entries  1 through  20 
determine  the  whole  matrix;  and  thus  the  storage  is  minimum. 

Given  an  element  (I,J)  of  the  persymmetric  matrix,  the  following 
algorithm  finds  the  location  K(I,J)  of  that  element  in  a vector  form. 


Let 

I 

In 

n,N/2} 

J 

in 

{I,(N+D  - 1} 

K 

in 

{l,N(N+2)/4} 

K(I,J) 

= { 

I-l 

I N - 2(il-l)}  + (J  - I)  + 1 
£=1 

Let 

I 

in 

{2,N/2} 

J 

in 

{1,(I-1)} 

K(I,J) 

= 

K(J,I) 

Let 

I 

in 

{2,N/2> 

J 

in 

{(N  + 2)  - I,N} 

K(I,J) 

= 

K[(N  + 1)  - J,(N  + 1)  - I] 

Using  the  Amdahl  computer  to  invert  the  convolution  matrix.  A,  the 
subprogram  NEWLOC.FTN  was  used  to  store  and  retrieve  the  entries  of  A”"^ 
in  a vector  form  to  and  from  the  department  computer,  POP  11/34.  The 

algorithm  to  restore  the  estimated  spectrum  is  presented  in  Appendix  A, 
Part  A-6. 
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In  the  case  of  symmetric  storage,  N-(N  + 1 )/2  memory  locations  were 
needed;  but  in  the  case  of  persymmetric  storage,  only  N*(N  + 2)/4  memory 
locations  were  needed.  For  the  case  N = 250,  the  necessary  storage  loca- 
tions are 

N*(N  + l)/2  = 31,375  in  the  symmetric  case 

N-(N  + 2)/4  = 15,750  in  the  persymmetric  case. 

Since,  in  double  precision,  the  highest  location  index  is  about  8,000  in 
the  POP  11/34,  the  preceding  algorithm  was  modified  to  accommodate  the 
department  computer.  The  Fortran  program  SY.FTN  is  an  implementation  of 
this  modified  version  and  requires  subprograms  SKIP.FTN,  KIJ.FTN,  and 
BIGEST.FTN.  The  running  time  of  1 hour  and  43  minutes  was  relatively 
very  good  compared  to  the  previous  cases. 

3.7  Comparison  of  the  Results  in  the  Cases  of  (50  x 501 
and  (250  x 250)  Convolution  Matrix  ^ 

Using  a sine  function  to  generate  the  50  x 50  convolution  matrix 
and  the  observed  spectrum  with  one  or  two  peaks  at  exact  or  slightly  dif- 
ferent frequency  locations,  the  results  were  as  expected  for  the  estimated 
spectrum  (i.e.,  sine  function  with  one  or  two  peaks). 

But  in  the  case  of  a (250  x 250)  convolution  matrix  generated  with  a 
sine  function  and  an  observed  spectrum  obtained  with  a Fast  Fourier  Trans- 
form, the  result  was  highly  noisy.  This  is  believed  due  to  round-off 
errors  in  the  inversion  of  the  large,  not-well -conditioned  matrix. 

3.8  Recovery  of  a Bandlimited  Function  Using  a Convolution 
Function  Obtained  from  the  Summation  of  Sine  Functions 

Figure  2.5  in  Chapter  II  suggests  that  there  is  a collective  effect 
among  the  regions  of  the  sampled  spectrum  (i.e.,  a given  region,  i,  is  not 
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independent  but  receives  a contribution  from  the  other  regions  of  the 
ensemble) . 

Therefore,  it  ought  to  be  noticed  that  the  convolution  matrix  should 
not  be  built  using  only  a single  sine  function  but,  rather,  a combination 
of  sine  functions  to  take  care  of  the  higher  order  spectra. 

Hence,  the  Fourier  Transform  of 

{100  zeroes  + sin(-^)  + 100  zeroes} 

I 

50  points 

is  not  a sine  function,  but  a combination  of  sine  functions  of  orders 

0,1,2,...  ,0°. 

A program,  SSA.FTN,  was  used  to  find  a sine  function  in  a given 
region,  i,  and  the  contribution  from  the  other  regions,  say  24  higher 
spectral  orders  as  in  Figure  2.6.  Following  the  usual  practice,  the  mean 
was  removed  from  the  data  before  processing.  The  Fourier  Transform  will 
yield  a zero  D.C.  level. 

With  the  assumption  that  the  estimated  spectrum  has  a Dirac  delta 
function  at  any  frequency,  the  observed  response  should  have  a D.C.  of 
zero  and  the  shape  of  a sine  function  with  a peak  value  at  the  same  fre- 
quency location  as  that  of  the  delta  function.  This  is  shown  in  the 
matrix  formulation  on  the  following  page. 

It  is  clear  that  moving  the  entry  1 in  the  vector  on  the  left  hand 
side  will  move  the  peak  of  the  sine  function  on  the  right  hand  side  of  the 
equal ity. 

The  determinant  of  such  a matrix  is  very  small,  the  product  of  the 
matrix  and  its  inverse  had  about  30  percent  error  from  the  identity  matrix. 

A . A"^  - I = ±kl, , k 0.30  . 
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The  fact  that  the  determinant  is  very  small  is  cured  numerically  by 
multiplying  A by  a constant  a to  get 

otA(n.n) 

which  means  that  the  determinant  is  multiplied  by  Using  program 
NCM.FTN  that  computed  a combination  of  sine  functions  with  30  spectral 
orders  250  points  away  from  one  another,  the  convolution  oA(n,n)  was 
built  and  inverted  with  the  university  computer  using  program  DCD.FTN 

Multiplying  the  result  by  a,  one  gets  Due  to  the  lack  of  funds  for 

the  university  computer,  it  was  decided  to  use  another  approach  to  solve 
the  problem. 
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3.9  Matrix  Inversion  Using  the  Toeplitz  Recursive  Algorithm 
3.9.1  Scheme 

The  convolution  matrix.  A,  beinq  Toeplitz  persymmetric,  one  may  take 
advantage  of  the  recursion  algorithm  for  the  solution  of  simultaneous 
equations  as  described  in  Part  7 of  Appendix  A.  Then  when 
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To  find  one  solves  for  the  first  column  of  A"^  such  as 
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a(ij) 

nxn 
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The  procedure  was  then  repeated  for  the  2nd  column,  3rd  column,  etc. 

3,9.2  Number  of  Operations  Involved  in  the  Solution  of  a Linear  System 
of  Equations  Using  the  Toeplitz  Algorithm  ~ 


n = 1 


n = 2,3,N 


calling 


1 

(n-1 ) 
1 

3n 

4(n-l) 

1 

2 

less  2N 
2(N-1) 
m: 
d: 
a: 
s: 
Op: 


multi  pi i cation 
division 

multiplications  for  'a' 

multiplication  for  'a' 

multiplications  for  '6,f,Y' 

additions  for  'a, 6,f,Y' 

subtraction  for  'q' 

divisions  for  'k  and  q' 

multiplications  for  last  6 and  y 

additions  for  last  3 and  y 

multiplication 

division 

addition 

subtraction 

the  number  of  operations. 
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N N N 

Op  = l(m)  + 1(a)  + I 4n(m)  + I 2(d)  + I 4(n-l)(a) 

n=2  n=2  n=2 

N 

+ I l(s)  - 2N(m)  - 2(N-l)(a) 
n=2 

Op  = l(m)  + 1(a)  + [2N(N+D  - 4](m)  + 2(N-l)(d)  + 2N(N-l)(a) 

+ (N-l)(s)  - 2N(m)  - 2(N-l)(a) 

reduced  to 

Op  = 2N^(m)  + 2N(d)  + 2N^(a  and  s) 
for  N = 250 

Op  = 124,997  multiplications  + 498  divisions  + 124,003  additions 
+ 249  subtractions. 

Subsequently,  to  obtain  the  inverse  of  the  convolution  matrix  with 
the  above  described  algorithm,  the  number  of  operations,  NOp,  is  N times 
bigger  than  Op  above.  And  to  find  the  deconvolved  spectrum,  NOp  is 
increased  by 

N[N(m)  + (N-l)(a)] 

Therefore  the  total  number  of  operations  needed  was 

^^otal  = 3(N^-l)(m)  + 2(N-l)(d) 

+ (3N^  - 4N  + l)(a)  + (N-l)(s)  . 

3.10  Sensitivity  of  the  Direct  Deconvolution  Technique 
To  investigate  the  sensitivity  of  the  direct  deconvolution  method, 
consider  the  case  where  the  basic  resolution  vector  is  obtained  from  the 
Fast  Fourier  Transform  of  a rectangular  window. 
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12  points 


18  points  18  points 


which  results  in  a sine  function. 

It  was  noticed  that  when  the  resolution  vector  and  the  observed  spec- 
trum are  identical,  the  result  is  good,  but  even  for  a slight  change,  say 
changing  the  zeroes  of  the  sine  function  from  0.000000  to  0.000002,  the 
result  is  not  good.  Thus,  the  presence  of  even  a minute  amount  of  the 
random  noise  destroys  the  inversion.  Instead  of  producing  the  desired 
delta  function  (identity  matrix),  the  original  sine  function  reappears. 

The  matrix  methods  for  direct  deconvolution  are  of  limited  value,  becoming 
increasingly  sensitive  to  noise  as  the  matrix  size  increases. 


CHAPTER  IV 

DECONVOLUTION  WITH  THE  "NIBBLING  TECHNIQUE" 

4,1  Development 

A second  deconvolution  method  is  studied  in  this  chapter.  Given  a 
spectrum  F(w)  convolved  with  a sine  resolution  function  R(w)  shown  in 
Figure  4.1,  the  Nibbling  technique  is  a cross  correlation  technique 
simplified  (i.e.,  minimization  or  elimination  of  the  multiplications 
involved  in  the  cross  correlation  operation). 

We  know  that  the  cross  correlation,  of  the  spectrum  F(w)  with 
the  resolution  function  R(w)  will  have  a maximum  value  where  the  two 
functions  are  highly  correlated.  Thus,  one  need  only  search  for  the  peaks 
of  F(w),  and  to  construct  the  true  spectrum,  Fj(w),  one  takes  a "bite" 
from  the  observed  spectrum  at  the  location  of  the  maximum  of  the  cross 
correlation,  ci)pp.  The  shape  of  the  bite  is  the  same  as  that  of  R(w).  The 
procedure  is  then  repeated  until  F(w)  has  been  "nibbled"  away. 

Because  of  the  multiplications  involved  in  the  cross  correlation 
operation,  it  will  undoubtedly  take  a longer  time  to  perform  than  the 
nibbling  technique  which  involves  just  additions. 

Indeed,  the  peak  of  the  observed  spectrum  corresponds  to  the  maximum 
of  the  cross  correlation  function  <()pp(w).  After  a reasonable  number  of 
iterations,  the  true  spectrum  is  recovered. 

The  positive  side  of  this  nibbling  technique  is  that  it  has  some 
virtues,  among  which  are  the  following: 
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Figure  4.1  Fourier  spectrum  F(w)  and  Resolution  function  R(w). 
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(1)  It  is  self-consistent. 

(2)  The  nibbling  of  a random  noise  is  close  to  zero,  as  expected 
from  the  cross  correlation  of  a sine  resolution  function  with 
a random  noise. 


R(w)  * Noise  = 0 

When  the  number  of  iterations  is  too  large,  instead  of  recovering 
a clean  spectrum,  a noise  contaminated  spectrum  will  result  as  depicted 
in  Figure  4.2. 

However,  the  method  also  has  a negative  side.  It  is  important  to 
indicate  that  this  deconvolution  method,  namely  the  nibbling  technique, 
can  be  used  only  in  spectra  involving  peaks;  and  this  is  almost  always 
the  case.  For  a filtered,  flat  spectrum  this  method  is  hopeless  and 
should  not  be  considered. 

4.2  Application  of  the  Nibbling  Technique 

4.2.1  Resolution  Function 

The  resolution  function  was  obtained  by  taking  the  Fourier  Transform 
of  a rectangular  window  with  80  zeroes  added  on  each  end  of  the  record. 

A correction  was  also  done  to  obtain  the  right  phase.  This  is  performed 
using  program  REC.FTN  with  subprograms  FOURG.FTN  and  COR.FTN. 

The  first  100  points  of  R(f)  constitute  the  resolution  function  as 
depicted  in  Figure  4.3, 

4.2.2  Observed  Spectrum 

The  observed  spectrum  is  the  Fourier  Transform 

F { 40  zeroes  + cos  [ ] + 40  zeroes  } 


20  points 


Clean  spectrum 
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The  spectrum  is  obtained  using  program  CP.FTN  with  subprogram  FOURG.FTN; 
the  result  is  also  rotated  to  the  right  phase  angle  using  program  COR.FTN 
and  is  depicted  in  Figure  4.4, 


4.3  Procedure 

Any  given  function,  f(t),  can  be  decomposed  into  an  even  function, 
fg(t),  and  an  odd  function,  f^Ct),  according  to 


f^(t)  = i(t)  * n-t) 
fjt)  = m - f(-t) 


Knowing  that  the  Fourier  Transform  of  a real  even  function  is  also 
a real  and  even  function 

F.T. 

and  the  Fourier  Transform  of  an  odd  real  function  is  an  imaginary  and  odd 
function, 

F.T. 

fo(t)  ► F^(w) 

then,  having  Fg(w)  and  F^(w),  these  spectra  are  nibbled  symmetrical ly  as 
follows. 

First,  the  maximum  value  of  Fg(w)  is  located  using  the  subprogram 
BLOC.FTN,  Then  translate  the  resolution  function,  R.j , along  the  frequency 
axis  to  make  the  maxima  of  Fg(w)  and  R-|(w)  coincide.  A second  resolution 
function,  R2(w),  is  located  at  a symmetrical  point  with  respect  to  the 
Nyquist  frequency,  Wj^^. 

At  this  point,  a new  resolution  function,  R(w),  was  created  using 
R-|  (w)  and  R2(w)  such  as 
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Figure  4.4  Observed  spectrum. 
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R(w)  = R^(w)  + R2(w) 

that  was  normal izGd  and  ussd  to  nibble  away  the  observed  spectrum  and 
recover  the  estimated  spectrum  Fg(w). 

Turning  to  the  odd  spectrum,  F^(w),  the  same  process  was  applied,  but 
in  this  case  instead  of  adding  the  resolution  functions,  R.j  (w)  and  R2(w), 
to  obtain  the  new  resolution  function,  R(w),  they  were  substracted  from 
each  other,  and  again  the  odd  spectrum  was  estimated  F^^(w). 

As  remarked  above,  Fg(w)  will  be  the  real  part  of  the  spectrum  F(w), 
and  Fq(w)  will  be  its  imaginary  part  and  both  were  used  to  find  the  magni- 
tude and  phase  of  the  estimated  spectrum  F(w)  as  follows: 

I F(w)  I = (F|(w)  + F“(w))'/^ 

Fo(w)  1 

KWJ  • 

J 

Or  else,  having  at  hand  a time  series,  f(t),  one  adds  zeros  on  each  end 
of  the  record  and  takes  the  Fourier  Transform  to  obtain  the  observed  spec- 
trum; then  locates  the  maximum  of  the  spectrum  and  nibbles  away  this 
observed  spectrum  until  the  remainder  is  small  enough.  This  point  is  dis- 
cussed in  the  following  section. 

4.4  Results 

To  begin  with,  four  things  are  defined. 

fraction:  is  that  part  of  the  resolution  function  that  is 

taken  away  from  the  observed  spectrum,  or  the 
"nibble"; 

e:  represents  the  largest  value  of  the  remainder; 

a:  standard  deviation  of  the  data  in  the  remainder; 


^ = tan 
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M.  thG  number  of  iterations  or  the  number  of  times  the 
resolution  function  was  used. 

Clearly  e,  a,  and  M are  the  parameters  that  will  terminate  the  nibbling 
process. 

Using  program  IBN.FTN  with  a resolution  function  close  to  the  observed 
spectrum  but  not  the  same,  for  different  values  of  the  convergence  criterion 
e,  say  e - 0.5,  0.35,  0.25,  0.15,  0.1,  0.075,  0.05,  the  number  of  itera- 
tions increases  from  69  to  558  for  e going  from  0.5  to  0.05,  respectively. 

However,  for  values  of  e = 0.35,  0.25,  ...  a leakage  in  the  form  of 
side  peaks  appears  in  the  deconvolved  spectrum  and  increases  slightly  with 
decreasing  e. 

The  same  process  was  repeated  using  a resolution  function  identical 
to  the  observed  spectrum  for  values  of  e starting  from  0.5  down  to  0.05, 
and  the  results  were  perfect— no  leakage  at  all.  And  the  number  of  itera- 
tions, for  the  same  value  of  e is  smaller  than  that  in  the  previous  case. 

As  an  example,  for  e = 0.05,  the  number  of  iterations  was  299  in  the 
present  case,  while  it  was  558  in  the  previous  case. 

When  the  resolution  function  and  the  observed  spectrum  have  the  same 
shape  but  not  very  close,  the  number  of  iterations  goes  up  and  there  is  a 
presence  of  leakage  in  the  deconvolved  spectrum.  For  e = 0.05,  M in  this 
case  goes  up  to  6,939. 

But,  when  another  criterion  is  used,  namely  the  standard  deviation 
of  the  data  in  the  remainder,  a,  to  end  the  process,  the  results  for  a = .15 
are  excellent  for  a cosine  wave  or  a sine  wave;  and  using  either  one  of 
the  nibbling  procedures  (i.e.,  symmetrical  or  not  symmetrical  nibbling). 

The  programs  used  were 

IBNN.FTN  for  the  symmetrical  nibbling 

NNBI.FTN  for  the  nonsymmetrical  nibbling. 
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In  the  case  of  widely  spaced  peaks  in  the  observed  Fourier  spectra 
(one  or  more  widths  of  the  main  lobe  of  the  resolution  function),  the 
nibbling  technique  of  deconvolution  gives  an  excellent  result,  namely 
delta  functions  at  the  precise  locations  of  the  peaks  in  the  observed 
Fourier  spectrum.  Thus  the  finite  width  resolution  function  convoluted 
with  the  true  Fourier  spectrum  is  replaced  by  the  true  Fourier  spectrum. 

When  observed  functions  were  studied  which  had  two  close  peaks,  the  results 
were  not  as  satisfactory.  As  the  peaks  came  closer  and  closer  together, 
the  deconvoluted  spectra  presented  satellite  peaks.  Instead  of  two  clean 
delta  functions,  we  obtained  four  peaks,  none  at  exactly  the  right  location. 
For  a separation  of  less  than  one  width  of  the  main  lobe  of  the  resolution 
function,  spurious  peaks  in  the  deconvoluted  Fourier  spectrum  appeared  and 
this  became  increasingly  worse  as  the  centroids  of  the  peaks  came  closer 
and  closer.  In  a further  test,  instead  of  searching  for  the  single  maximum 
of  the  observed  spectrum  and  nibbling  at  that  point,  the  largest  sum  of  3, 

4,  and  5 points  was  searched  for  first,  and  the  nibbling  performed.  This 
procedure  reduced  the  number  of  spurious  peaks  and  improved  the  deconvoluted 
spectrum. 

The  methods  above  were  equivalent  to  performing  a cross  correlation  of 
the  resolution  function  with  the  observed  Fourier  spectrum  in  order  to  find 
the  point  of  maximum  correlation  and  then  subtracting  a nibble  centered  at 
that  point.  Only  a single  point  of  the  resolution  function  was  used  for  the 
cross-correlation  however.  It  is  felt  that  significantly  better  results  may 
be  obtained  if  a larger  portion  of  the  resolution  function  is  used  in  the 
cross-correlation,  perhaps  all  of  the  points  in  the  central  lobe,  or  even 
including  some  of  the  side  lobes.  (Note  that  the  mean  of  the  Fourier  spectrum 
should  be  removed  before  each  cross-correlation).  This  method  should  permit 
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deconvolution  of  closely  spaced  peaks,  but  would  clearly  consume  much 
more  computer  time  due  to  the  number  of  multiplications  required. 


CHAPTER  V 

EXTRAPOLATION  OF  A BANDLIMITED  FUNCTION 
WITH  AN  ITERATIVE  TECHNIQUE 

5,1  Introduction 

It  is  well  known  that  a bandlimited  function,  f(t),  is  an  entire 
function  and  can  be  expanded  in  a Taylor  series  as 

f(t)  = ! f (5.1) 

n=-oo 

that  converges  for  any  t [8], 


-T  0 +T 

The  derivatives,  f^'^^(o),  are  determined  from  f(t)  in  a finite  segment 
about  the  origin  called  f^(t).  In  principle,  this  provides  a method  to 
extrapolate  f.^(t).  But,  first  of  all,  there  is  difficulty  in  evaluating 
the  derivatives  of  f(t)  about  the  origin  in  the  presence  of  noise,  and 
secondly  a truncation  is  also  introduced  in  the  computation  f(t)  from 
Equation  5.1  above.  Errors  in  the  evaluation  of  f(t)  are  fatal  for  the 
determination  of  the  spectrum  F(u3). 

As  indicated  in  Chapter  II,  an  iterative  method  introduced  in  1975 
by  A.  Papoulis  could  be  used  to  extrapolate  a given  function,  f(t),  based 
on  the  Fourier  Transform. 

The  scheme  consists  of  the  following: 

(1)  Determine  the  Fourier  Transform  F^(co)  of  f.^(t)  in  a first  itera- 
tion noted  by  index  1, 
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(2)  Find  F^(w)  = F^(to) .p^(to) 

(3)  Next  the  inverse  Fourier  Transform  of  F^(o))  is  taken 

fl(t)  = r^(F^(o))) 

to  form 

gi(t)  = f^(t)  + [f(t)  - f^(t)]  p^Ct) 


g-|(t) 


f^(t)  , |t|  < T 

f-|(t)  , |t|  > T 


(4)  Then  the  direct  Fourier  Transform  is  taken 
G^(w)  = F(g^(t)) 

At  iteration  n 


Fn(“)  = 

(5)  Then 

f„(t)  = F-'  { F^(o.)  ) 
and 

9„(t)  = f„(t)  + [f(t)  - p^(t)  . 

The  extrapolated  spectrum  is  given  by 
GnM  = F {g^(t)> 

After  a number  of  iterations  F^(o))  will  approach  the  sought  after 
extrapolated  spectrum  F(o)) 

lim  F (o))  = F(o3) 

n-Mx. 
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5.2  Examination  of  the  Iterative  Technique 

To  start,  a sampled  time  series  was  obtained  as  follows: 

f^(t)  = 60  zeroes  + cos(^^g~^  + 60  zeroes 

t 

30  data  points 

Then  the  Discrete  Fourier  Transform  was  taken  to  obtain  F^(cj)  as  depicted 
in  Figure  5.1. 

To  multiply  the  spectrum  F (co)  by  the  truncation  window  p (w)  is 
equivalent  to  zeroing  the  spectrum  from  (oj^  - (z/2))  to  + (z/2))  as 
depicted  in  Figure  5.2.  Taking  the  inverse  Fourier  transform  of  F^(o)) 
should  give  f-|(t)  as  shown  in  Figure  5.3.  And  then  taking  once  again  the 
Fourier  Transform  of  g^(t)  to  get  &,(cj)  the  procedure  may  be  repeated  a 
number  of  times  to  recover  the  extrapolated  spectrum. 

This  scheme  has  been  performed  using  program  NTE.FTN  with  subprogram 
ITERA.FTN  for  different  lengths  of  the  section  of  zeroes  above  and  below 
the  Nyquist  frequency  and  25  iterations.  It  is  noticed  that  the  cosine 
wave  shape  is  recovered  but  the  extrapolated  section  has  a relatively 
small  amplitude. 

A comparison  of  the  extrapolation  with  different  numbers  of  iterations 
(e.g.,  25  and  101)  was  performed  and  the  results  were  the  same.  The  organi- 
gram is  given  in  Figure  B.l  of  Appendix  B. 

Improvement  of  resolution  in  the  frequency  domain  is  equivalent  to 
extrapolation  in  the  time  domain.  Our  results  indicate  that  the  iterative 
method  does  not  always  converge  in  amplitude  to  the  desired  function  even 
in  the  absence  of  additive  noise.  The  presence  of  such  noise  would  greatly 
increase  the  uncertainties  in  the  time  extrapolation. 


60  zeroes  I \ / \ / I 60  zeroes 
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Figure  5.3  First  iteration  (b). 


CHAPTER  VI 

EXTRAPOLATION  OF  A BANDLIMITED  FUNCTION  WITH  USE  OF 
THE  PROLATE  SPHEROIDAL  WAVE  FUNCTIONS 


6.1  Motivation 

It  is  possible  to  improve  the  frequency  resolution  of  the  Fourier 
Transform  (FT)  of  a timelimited  function  by  extrapolation  of  that  func- 
tion followed  by  Fourier  Transformation  of  the  result.  It  has  been 
suggested  that  this  is  possible  through  use  of  Prolate  Spheroidal  Wave 
Functions  which  have  uniquely  appropriate  properties:  they  are  orthogonal 

on  two  intervals. 


and 


where 


-00  < t < +«° 


-T  < t < +T 


r 

♦„(t)  *„(t)  dt  = x„ 

(fi.l) 

-T 

/.+  oa 

-oo 

♦„(t)  *„(t)  dt  = 

(6.2) 

<fp(t)  = prolate  spheroidal  eigenfunction 

= eigenvalue  corresponding  to  4>p(t) 

= Kronecker  index 
mn 


Thus  if  f(t)  is  known  on  the  interval  -x  < t < t and  is  Q bandlimited, 
then  the  extrapolated  version  of  f(t)  may  be  written: 


n=0 


a * (t) 
n ^n'  ' 
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with 


a 


n 


♦n't)  • dt 


— 00 


However  because  of  Equation  6.1 

r +T 


a 


n A 


*„(t)  dt 


f + T 


f(t)  4)^(t)  dt 


n 


-T 


and  so 


A summary  of  the  properties  of  these  functions  is  contained  in  Appendix  C. 
The  remainder  of  this  chapter  is  devoted  to  the  (difficult)  problems 
involved  in  calculating  and  handling  these  functions. 

The  functions  4>p(t)  above  are  purposely  normalized  versions  of  the 
angular  prolate  spheroidal  functions,  solutions  to  the  equation: 


which  is  obtained  when  expressing  the  wave  equation  in  a prolate  spheroidal 

coordinate  system.  The  X„  are  eigenvalues  chosen  so  that  the  S (c,t) 

n ' mn 


which  (as  we  will  see)  is  related  to  the  frequency  bandwidth  of  the  func- 
tion to  be  extrapolated. 


In  this  section  the  eigenvalues  X„(c)  and  the  eigenfunctions  S (c,t) 

mn  mn 

in  the  region  -1  < t < +1  are  determined. 

6.2.1  Expansion  in  Terms  of  a Series  of  Pn(t) 


^mn^^’^^  remain  finite  at  the  singular  points  t = ±1.  The  c is  a parameter 


6.2  Eigenfunctions  and  Eigenvalues 


To  determine  the  eigenvalues  of  the  prolate  angular  spheroidal  wave 
function  equation  for  the  case  m = 0 and  c = 0 (Legendre's  equation) 
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(6.3) 


whose  solution  is 


with  eigenvalues 


and  with  the  boundary  condition  " ^n^°^  ’ 


To  determine  the  number  of  mesh  points  necessary  for  an  accurate 
solution  for  the  eigenvalues  in  the  case  of  c not  zero,  one  varies  the 
number  of  mesh  points  in  the  case  c = 0 and  solves  for  the  eigenvalues 
which  are  known  to  be  n(n  +1).  The  number,  N,  of  mesh  points  that 
will  give  a result  very  close  to  n(n  + 1)  will  be  chosen  to  find  the 
eigenvalue  in  the  case  c not  zero  for  the  same  order  n. 

To  solve  the  angular  prolate  spheroidal  wave  equation  numerically 
for  the  eigenvalues 


for  m = 0 


(6.4) 


Let 


u,  = S 


U2  = = u 


Equation  6.4  becomes 


u 


Hence 


u 


This  is  a system  of  two  differential  equations  with  boundary  conditions 
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u,,„(0)  = P„(o) 
and 

“l,n<°>  = “ 

This  system  of  differential  equations  was  solved  using  a Runge-Kutta 
of  order  four  for  4 points  and  then  using  a Adams-Moulton  technique  for 
the  remaining  points  to  obtain  an  accurate  solution.  The  programs  used 
are  NUMEIG.FTN,  INITI.FTN,  R.FTN,  and  ADAM.FTN  and  are  available  in  the 
Department  of  Nuclear  Engineering  Library  at  the  University  of  Florida. 
6.2.2  Asymptotic  Expansion 

For  large  values  of  the  space  bandwidth  parameter  c and  higher  order 
n,  the  following  asymptotic  expansion  as  developed  by  Meixner  [9]  using 
the  parabolic  cylinder  functions  and  their  recursion  formulas  was  used  to 
obtain  the  eigenvalues. 

X„,_(c)  = (2n  - 2m  + 1)  c 

C-K» 

A more  accurate  limit  for  m = 0 is 

Tim  = (2n  + l)c  - (2n^  + 2n  + 3)2'^ 

C-KX) 

- (2n  + D(n2  + n + 3)2" V 

- [BCn**  + 2n'  + 7n  + 3)]2'®  c"^ 

- [66n®  + 1650**  + 962n^  + 1278n^  + 1321n  + 453]2"^°c"^ 

- [252n®  + 756n®  + 5885n“  + 10510n^  + 18478n^ 

+ 13349n  + 4425]2"^%"“ 
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- [527(2n  + 1)’  + 61529(2n  + 1)®  + 1043961  (2n  + 1)^ 
+ 2241599(2n  + 1)]  2"^°  c’^ 


The  program  used  is  EIGEN. FTN. 


6.2.3  Expansion  in  Legendre  Polynomials 

At  this  point,  the  eigenvalues  are  found.  To  determine  the  solution 
Son(c.t)  one  needs  to  find  the  expansion  coefficients  of  the  infinite 
series,  d®*^  (c)(see  Appendix  C),  for  the  case  m = 0. 

The  recurrence  formula  for  these  expansion  coefficients  is  given  by 
Equation  C.7  where  m is  zero. 


d°+2(c)  + [r(r  + 1 ) 


(c)  + 


2r(r  + 1)  - 1 
r^r  - lj(2r  + 3) 


c^] 


• + r(r  - l)c^  jOn  , X 

W - 3)(2r  - 1)  ^2^^^)  " ° 


and 


= L d°"(c)  P ft) 


on 


r=0 


from  Equation  6.5  the  d^_j_2  coeffiC-ients  are  given  by 
" |r'Vl)|^^^  2)c2  -r(r+l) 


(6.5) 

(6.6) 


2r(r  + 1)  - 1 
Cir  - 1)(2r  + 3) 


c^) 


(r  - 1 )r  c^ 

3)(2r  - 1) 


(6.7) 


For  order  n even  the  Equation  6.6  becomes 

S (c,t)  - I d (c)  P (t)  , r even  > 0 (6.8) 

r=0,2,4,...  ^ 


Let  r = 0 in  Equation  6.7  to  get 


^ dj(c)  - [X„(c)  - si  ] dS(c)  = 0 

d^(c)  = a^(c)  dS(c)  -f 

Let  r = 2 in  Equation  6.7 

“4('>  ° T§^  ® - TT) 

= d3(c)  d5(c) 

where 

a;(c)  = ilfr  [(X„(c)  - 6 - 1^)  a^(c)  - | ] . 

Thus  writing 

dj(c)  = a|;(c)  dj(c) 

with 

a3(c)  = 1 

where  the  new  expansion  coefficients  aJJ{c)  were  obtained  through  the 
same  recurrence  formula  as  the  djl(c)  coefficients.  And  replacing  the 
djl(c)  coefficients  by  aj(c)  dJCc)  in  Equation  6.8  yielded 

Sp(c.t)  = d^  a|:(c)  . P^(t)  . (6.9) 

Therefore,  for  any  given  pair  of  values  c and  t,  the  summation,  I, 
above  can  be  calculated  as  a function  of  the  eigenvalue  X|^(c)  corresponding 
to  the  eigenfunction  S^(c,t). 

For  an  odd  order  n.  Equation  6.6  becomes 
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Let  r = 1 in  Equation  6.7 

^ 6c2  ^ 

= b^Cc)  d!^(c) 

with 


Let  r = 3 in  Equation  6.7 

<^sM  =5^[(>!n(c)  - 12  -§c^)  d5(c)  ] 
or 

d"(c)  = b^(c)  . d"(c) 

with 

'’s*')  = ^ - 12  - i ^ 3 ■ 

Hence  writing 

d"(c)  = b;i(c)  • d!|(c) 

with 

b^Cc)  = 1 

and  the  new  expansion  coefficients,  bjl(c),  are  obtained  from  Equation  6.7. 
Equation  6.10  becomes 

SpCc.t)  = d!|'(c)  • b[l(c)  P^(t)  (6.11) 

Again,  the  summation,  I,  in  the  above  expression  can  be  calculated  for  any 
given  pair  (c,t)  as  a function  of  the  eigenvalue. 

In  both  cases  (n  even  and  n odd)  one  must  find  the  set  of  coefficients 


for  n even 
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and 


with 


{b^j  b^j  bgj,,,,b|^^}  tor  n odd 


Let  us  note  that  the  recurrence  relationship  (Equation  6.7)  consti- 
tutes a homogeneous  difference  equation  of  the  second  order  and  corresponds 
to  a differential  equation  of  the  2nd  order  and  therefore  has  two  nontrivial 
independent  solutions. 

Let  us  find  the  limit  when  r goes  to  infinity  of  Equation  C.7. 


r c 


.mn 


4r2  r+2 


+ [r' 


2r^ 


c"]  d™  + d™, 

4r2  J 4^2 


= 0 


or 


c .mn 
4 °r+2 


+ r^  d"’" 


X c jmn  _ „ , . 


(6.12) 


(6.13) 


We  have  two  cases  at  hand;  either  the  coefficients  dj!*^  constitute  a 


convergent  expansion  or  a divergent  one, 


jmn  .mn  , c^  .mn  ^ 
''r+2  * '•r  4F5-  <<r-2  ° 


(6.14) 


(i)  for  a convergent  expansion  we  have 


jmn  .mn  . .mn 
dr+2  < ^r  < '^r-2 


^nin 

^p+2  ^9  ^ 

and  2.  9oes  to  zero  faster  than  and 

r^  y.2  ^ 

And  consequently  the  first  term  in  Equation  6.14  is  neglected 
to  get 


(6.15) 


(ii)  for  a divergent  expansion  we  have 


.mn 

> d'""  > 
r 

,mn 

^-2 

.mn 

''r+2 

> d™  > 

jmn 

'^r-2 

r2 

r^ 

r^ 

case 

the  third 

term  in 

4r2  ^ 

imn  . ,mn 
'r+2  * ='r 

= 0 

(6.16) 


Hence 


.mn 


Equation  6. 15  — lim 


r->oo  d 


mn 

r-2 


4^  = 0 (convergent  expansion) 


the  [ratio  1 <1  , the  series  is  absolutely  convergent  [15]  (6.17) 

jmn 

4r  ^ 

“ ~5T  “ " (divergent  expansion) 


Equation  6.16 


1 im  ^ 


r-x»  d 


mn 

r-2 


Therefore  one  must  be  careful  in  the  determination  of  the  expansion 
coefficients. 

As  a matter  of  fact  when  trying  to  determine  the  expansion  coefficients 
3j„(c)  and  b^(c)  using  the  general  form  given  by  Equation  6.7  as  written  in 
program  SUBEXP.FTN  one  automatically  hits  the  divergent  series.  To  avoid 
this  problem  it  is  better  to  start  from  the  last  term  in  the  expansion 
and  continue  downwards  using  the  following  form  rather  than  Equation  6.7. 


r-2 


- I2r  - 3)(2r  - 
r(r  - 1)  c= 


11 


[(x  - r(r  + 1) 


2r(r  + 1)  - 1 
T2r  - I j(2r  +3) 


c^) 


In  this  form  one  must  input  a guess  for  d'^(c)  and  d"  „(c) 

r r+2'  ' 
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(r  + 2 = Nq  , r = Nq-2  for  n even,  and 

r + 2 = , r = N^-2  for  n odd). 

However  d^^^  or  dj^^  and  d^^-2  or  dj^^-2  should  be  small  to  avoid  an  overflow 

because  the  coefficients  will  be  larger  and  larger.  And  once  the  a*^(c)  or 

r 

bj^(c)  are  found  they  should  be  normalized  to  a^Cc)  = 1 and  b!|'(c)  = 1. 

In  this  manner  the  new  coefficients  are  finally  found: 

C 1 , » • • • » 3 

or 

r I5  bg,...,b|^_^  J 

To  recover  the  expansion  coefficients  d[I(c)  the  above  sets  of  coeffi- 
cients must  be  normalized  as  follows.  For  n even  we  have 

s„(c.t)  = r d"(c)  pjt) 

= dj(c)  1°  a"(c)  P (t) 

u r=0  ^ 

The  normalized  constraint  is 

Lim  S (c,t)  = S (c,o)  = P^(o) 
t->0  '>  n n 

Thus 


For  n odd  we  have  a similar  expression 

s„(c,t)  = d!;'(c)  r b”u)  p^(t) 

r=l 


where 
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But  in  this  case  (n  odd)  there  is  a small  difficulty  because  the 
Legendre  polynomials  at  t = 0 are  all  zeroes  and  the  coefficient  d!j'{c)  will 
be  undetermined  (i.e.,  d!j'(c)  = (0/0). 

To  find  d^(c)  one  way  is  to  use  Hospital's  theorem  to  evaluate  the 
limit  of  d!j'(c)  when  t goes  to  zero.  Hence 


and  we  use  the  property  of  the  Legendre  polynomials 


to  get 


Pp(o)  = r P^_-|(o)  , 


(o) 


r • 


(o) 


(6.20) 


Once  dQ(c)  and  d-j(c)  are  obtained  via  Equation  6.18  or  Equation  6.20, 
respectively,  the  expansion  coefficients  d|](c)  are  completely  determined. 
The  above  scheme  was  applied  using  the  following  programs. 


PRO.FTN 

SUB.FTN 

LEGR.FTN 


as  main  program 

as  subprogram,  subroutine  to  find  the 
substitute  coefficients  ajl(c)  or  bjl(c) 
subprogram,  subroutine  to  find  the 
Legendre  polynomials. 
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Hence  the  parameter  c as  given  by  Equation  C.17  becomes 


(6.22) 


Finally,  given  a record  of  N data  points,  one  determines  the 
parameter  c from  Equation  6.22  and  finds  the  eigenvalues  and  eigenfunctions 
as  indicated  above  and  proceeds  to  expand  the  bandwidth  limited  function, 
f(t),  in  terms  of  the  prolate  spheroidal  wave  functions  <i)^(c,t). 

N 


where 


f(t)  = I a 4)  (c,t) 
n=0  " ^ 


®n  " 


knowi ng 


T/2 
-T/2 


f(t)  cf)^(c,t)  dt 


S^(c,t) 


with  S^(c,t)  the  angular  prolate  spheroidal  function  and 

+1 


4m  = 


Son  ''t 


■1 


Instead  of  calculating  separately  the  expansion  coefficients  a^,  it  is 


To  minimize  the  error,  however,  one  should  use  a compact  scheme. 
;ad  of  calculating  separatel 
preferable  to  use  the  following. 

f(t) 


N 

I 

n=0 


/ 1 

1 ^ 

'■  n 

. 

T/2 


X 


■T/2 


f(t)  [ — ^ S (c,t)  dt 
n 


s„(c,t) 

n 


f(t)  = 


N 

1 

n=0 


T/2 

■T/2 


f(t)  S^(c,t)  dt 


S^(c,t) 
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This  concludes  the  recovery  of  the  bandlimited  function  f(t)  in  the 
interval  -1  < t ^ +1  where  in  the  above  formulation  T = 2. 

6.2.4  Expansion  in  Terms  of  Spherical  Bessel  Functions 

Outside  the  interval  -1  ^ t <:  +1 , the  angular  prolate  spheroidal 
function  is  not  defined  in  terms  of  the  Legendre  polynomials. 


For  t in  [-1,  +1] 
N 

I 

r 


Sn(c.t)  = 1°’'  d|;(c)  P^(t) 

"l 


(6.23) 


and  the  radialprolate  spheroidal  function  of  the  first  kind 

,0) 


Rp  (c,t)  = k^(c)  S^(c,t) 


(6.24) 


and  for  t in  -1]  U [+1,  +®] 

d;:  • 1“’’  i'--"  dj(c)  z<')(ct) 


where 


zJi  ^ct)  = (2^^)  '^n+(l/2)^'^^^ 


and 


" ^2ct^^^^  '^n+(l/2)^‘^^^ 


consequently 


= spherical  Bessel  function  of  order  n 

d^(c)  j„(ct) 


(6.25) 
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Thus,  one  should  evaluate  Equation  6.24  for  t = 1 and  Equation  6.25  for 
t = 1 and  use  the  constraint  that 

R*”(c,n  = Rt’>(c,r) 


1 = t approaching  from  below 
= t approaching  from  above 

Equation  6.24  yields 

R^^^c.D  = a 

Equation  6.25  yields 

= B 


Let  c.  = N„6  N„  = | 

Therefore,  one  must  find  the  radial  prolate  spheroidal  function 

outside  the  interval  [-1,+1]  using  the  spherical  Bessel  functions  and 

multiply  the  result  by  the  normalizing  factor  N . 

n 

The  angular  prolate  function  is  then  found  using  the  joining  factor 
such  that 

S^(c,t)  = K^(c)  • R|^"*hc,t)  for  |t|  > 1 
The  spherical  Bessel  function  is  given  by 

Jn(z)  = f^(z)  . sin(z)  + (-l)"'*’^  ^"-(n+l)(^)  • 
where  the  f^'s  are  elementary  functions  following  the  recurrence  formula 


Vl^^)  ^n4.l(2) 


n+1 


(2n  + 1)  f^(z) 
z 


for  n = 0,  ±1,  ±2,..., 
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with 

fgCz)  = 1/Z 
f^(z)  = 1/z^ 

This  scheme  was  applied  using  programs: 

FXTM.FTN 

SUB.FTN 

LEGR.FTN 

SUMAP.FTN 

BESSEL. FTN 

SUMBES.FTN 

SUMD.FTN 

JOIFAC.FTN 

FACT. FTN 

BIGEST.FTN 

in  the  library  of  the  Department  of  Nuclear  Engineering  of  the  University 
of  Florida. 

6.2.5  Conclusion 

The  preceding  methods  were  applied  to  extrapolate  a bandwidth  limited 
function  in  the  form  of  a cosine  wave  in  the  interval  [-!,+!].  The 
result  was  as  depicted  by  Figure  6.1.  Inside  the  interval  [-1,+1]  the 
function  was  recovered  but  distorted;  on  the  other  hand,  outside  this 
interval  there  was  a predominance  of  the  high  order  prolate  functions  and 
the  shape  is  right  but  the  frequency  is  higher,  about  twice  as  large. 

We  concluded  that  the  methods  described  thus  far  were  not  adequate 
for  sufficiently  high  orders  n.  High  orders  are  necessary  for  extrapolation 
purposes,  and  the  expansions  of  SQ^(c,t)  in  sums  of  Legendre  polynomials 
or  parabolic  cylinder  functions  were  not  adequate  for  this  purpose. 

6.3  Frobenius  Expansion 

It  was  decided  to  try  a power  series  expansion  around  the  singular 
point  at  t = 1 . The  resultant  S^^(c,t)  would  certainly  converge  in 
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-1  ^ t < 3 and  enable  at  least  a two  fold  extrapolation.  The  boundary 
conditions  used  were: 

S"(c,l ) _ x-c^ 

S(c,l)  2 

and 

S^(c,0)  = 0 or  S(c,0)  = 0 
depending  on  the  order  n. 

6.3.1  Development  of  the  Method 

Rewriting  the  differential  equation, 

(1  - t^)  ^ S - 2t  ^ + (X  - c^t^)  S = 0 (6.26) 


let 


P(t)  - (1 


q(t)  = 


X - C^t^ 

(T  - F) 


the  equation  becomes 


0 + P(t)  ^ + qCt)  S = 0 [16] 


The  differential  equation  has  a singularity  at  t = ±1,  and  is  solved  using 
a Frobenius  expansion  about  t = +1 , 

s = I a (t  - 1)"*'’ 

n=0  " 


Z a (n  + v)(t  - 
n=0  " 

0 “ l a (n  + v)(n  +v-l)(t  - 
n=0  ^ 
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Bringing  these  expressions  into  Equation  6.26 

(1  - t^)  I (n  + v)(n  + V - 1)  a (t  - . 2t  I (n  + v) 

n=0  " n=0 


,n+v-l 


a (t  - 1)"^”-'  + (X-  I a (t  - 1)"*''=  0 

^ n=0  ^ 


(6.27) 


Let 


z = t - 1 


t = z + 1 

(1  - t^)  = (1  - t)(l  + t) 
= -z(z  + 2) 


Equation  6.27  becomes 

-z(z  + 2)  (n  + v)(n  + v - 1)  a - 2(z  + 1)  I a (n  + v) 

n=0  ^ n=0  " 


.n+v-1 


( X-  c^z^  - 2c^  z - c^)  I a = 0 

n=0  " 


Grouping  terms  of  the  same  power  of  z yields 


I [(n  + v)(n  + V - 1)  + 2(n  + v)  + (c^  - x)]  a z 


n+v 


n=0 


+ I [-2(n  + v)(n  + V - 1 ) - 2(n  + v)]  a z 
n=0  ^ 


n+v-1 


+ I [-2c^]  a„  X I [-c“]  a„  = 0 

n=0  ^ n=0  " 


(6.28) 


zn+v-  1 ^ coefficient  -2(n  + v)(n  + v-  1 +1)  a. 


.n+v 


has  a coefficient  -[2(n  + v)  + (c^-  X)  + (n  + v) 


(n  + V - 1 )]  a. 


80 


^n+vH-l 

has  a coefficient 

-2c 

^n+v+2 

has  a coefficient 

-c2 

For  Equation  6.28  to  be  satisfied,  each  coefficient  of  a power  of  z 
should  be  zero.  Thus  for  the  power  of  z,  say  (n+\H-2)  one  should  have 

-c^  a^  - 2c^  a^^.|  - [(n  + v + 2)(n  + v + 2 - 1 ) + 2(n  + v + 2) 

+ (c^  - x)]  - 2(n  + V + 3)"  a^^3  = 0 (6.29) 

Let  n = -3  with 

a_3  = a_2  = a_^  = 0 
Equation  6.29  yields 

-2(-3  + V + 3)2  a^  = 0 
v2  = 0 

V = 0 is  a double  root. 

Equation  6.29  becomes 

-2(n  + 3)"  - [(n  + 2)(n  + 1 ) + 2(n  + 2)  + (c“  - x)] 

-2c“  = 0 

Let  m=n  + 3 n = m-  3 to  get 

-2m2  - [(m  - 1 ) (m  - 2)  + 2(m  - 1 ) + (c^  - x]  a^_.j 


-2c2 


a^  , - c2  a„  - = 0 
m-2  m-3 


(6.30) 
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Let  k=m-lm=k+l 

-2(k  + 1)^  aj^^^-[k(k  - 1)  + 2k  + (c^  - x)]  aj^  - 2c^  a|^_^ 

- a,^_2  = 0 (6.31) 

Rewriting  this  equation  for  successive  values  of  k. 


0 

II 

to  get 

-2  a.|  - (c^  - x)  ag  = 0 

(6.32) 

k = 1 

to  get 

-2  X 4 X 82  - [2  + (c^  - x)] 

a-|  - 2c^  ^0  ~ ® 

(6.33) 

k = 2 

to  get 

-18  83  - [6  + (c^  - x)]  a2  - 

2c^  a-j  - c^  8q  = 0 

(6.34) 

Equation  6.32  yields  a^  = -(c^  - x)  3q  / 2 
Equation  6.33  yields  82  = -(2c^  + (2  + (c^  - x))  a-j)  / 8 
Equation  6.34  yields  a^  = -(c^  + 2c^  a^  + (6  + (c^  - x))  ^ 

and  Equation  6.31  yields 

[k(k  - 1)  + 2k  + (c^  - x)]  a,^  + 2c^  a|^_-,  + c^  a,^_2 


Therefore,  one  sets  8q  = 1 and  finds  the  other  coefficients 
ag  = 1 

a^  = (x  - c^)  / 2 

a2  = [((x  - c^)  - 2)  a^  - 2c^]  / 8 

- [(k(k  - 1)  + 2k  + (c^  - x))  aj^  + 2c^  a^_-^  + c^  a,^_2] 
^k+1  " ’ 


2(k  + 1)^ 
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Once  the  expansion  coefficients  are  computed,  the  angular  prolate 
spheroidal  function  of  order  n is  given  by 


S (c,t)  = I a.  (t  - 1)' 
" k=0 


or 


S (c,t)  = aQ  I ~ it  - 1)*^ 

" ^ k=0  ^0 


The  coefficient  aQ  is  obtained  from  the  normalization  constraint, 

Sj^(c,o)  = P^(o) 

where  Pp(o)  is  the  Legendre  polynomial  of  order  n evaluated  at  zero.  For 
even  order  n,  assuming  aQ  = 1 , the  constraint  becomes 


l a,  (-1)"  = P„(o) 


and 


k=0 


a«  = 


P„(0) 


° n If 

k=0 


for  odd  order  n,  Pp(o)  = 0,  aQ  will  be  determined  using  Hospital's  theorem 


S'(o)  = I a k(0  - 1)' 
" k=0  ^ 


and 


=0  ' 


I 

k=0 


Thus 
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S^(c,t)  = 


Pjo) 


I 

k=0 


Z aJt-i)'' 

k=0 


n even 


S (c,t)  = I a.  (t  - 1)^  , n odd 

I • k . (-1)^  k=0 

k=0 

This  Frobenius  expansion  about  t = 1 is  valid  until  the  next  singular 
point  is  reached,  namely  t = -1.  Therefore  the  radius  of  convergence,  R, 
of  the  solution  is  two. 


-3 


2 


3 


This  extrapolation  of  a bandlimited  function  can  be  performed  up  to 
t = + 3 as  follows 

f(t)  = I a^  <J>p(c,t) 
n 

with 

<!>n(c,t)  = S^(c,t) 


where  4'p(c,t)  is  the  prolate  spheroidal  function.  The  coefficient  a^  is 
obtained  using  the  orthogonality  property 


and 


n 


+1 

f(t)  (})^(c,t)  dt 

. -1 


f(t) 


n n 


f(t)  S^(c,t)  dt  ] S^(c,t) 


-1 
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where 


+1 

<})^(c,t)  dt 


+1 

S^(c,t)  dt 

. -1 


S^(c,t)  dt 

-1 

and  f(t)  becomes 

f(t)  S^(c,t)  dt 

f(t)  = I • S (c,t) 

n +1 

S^(c,t)  dt 

. -1 

6.3.2  Results 

The  above  described  scheme  was  applied  with  the  following  input 
functions  in  the  domain  -1  < t < +1  in  order  to  find  the  values  of 
c and, n. for  which  extrapolation  was  possible.  - 

(a)  f(t)  = $2(5, t) 

(b)  f(t)  = $3(5, t) 

(c)  f(t)  = S^(5,t)  + 82(5, t)  + 83(5, t) 

(d)  f(t)  = cosine  function  with  6 zeroes  in  [0,+l] 

(e)  f(t)  = cosine  function  with  1 zero  in  [0,+l] 

(f)  f(t)  = cosine  function  with  1 zero  in  [0,+l] 
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In  cases  (a),  (b),  and  (c)  the  result  is  excellent,  the  input  is 
recovered  perfectly. 

In  case  (d)  where  the  space  bandwidth  parameter  c = 5 and  using  7 
eigenvalues,  clearly  the  result  is  bad  as  expected  because  at  least  12 
eigenvalues  are  needed  to  recover  the  input. 

In  case  (e),  with  7 eigenvalues  and  c = 5,  the  result  was  very  good. 

In  case  (f),  which  is  the  same  as  case  (e)  but  where  c = 50,  the 
result  is  very  bad.  This  is  due  to  the  fact  that  the  expansion  coeffi- 
cients, aj^,  are  truncated  while  they  are  still  very  large. 

To  determine  where  to  stop  the  expansion,  for  the  following  values 
of  the  space-bandwidth  parameter  c and  using  the  asymptotic  expansion  to 
generate  the  eigenvalues  x^  with  a negligible  residual  it  is  found  that 
the  coefficients  aj^  are  still  large  and  for  large  values  of  c they  even 
give  an  overflow. 


c = 

10 

7 eigenvalues  used 

c = 

25 

7 eigenvalues  used 

c = 

50 

7 eigenvalues  used 

c = 

100 

7 eigenvalues  used 

and  k = 90 
max 

= 125 
= 175 

the  coefficients  aj^  diverge. 


It  was  concluded  that  this  method  requires  one  to  retain  a large 
number  of  terms  in  order  to  extrapolate,  and  that  the  eigenvalue  expansions 
used  to  calculate  the  Xp  were  inadequate  for  this  purpose.  It  was  decided 
to  examine  direct  numerical  integration  of  the  differential  equation. 


6.4  Direct  Numerical  Integration  of  the  Differential  Equation 
In  this  section  the  differential  equation  for  the  angular  prolate 
spheroidal  function  is  solved  numerically  using  a Runge-Kutta  and  explicit 
Adams-Moul ton  method. 
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6.4.1  Boundary  Conditions 

We  know  that  the  angular  prolate  spheroidal  function,  S^(c,t),  and 
its  derivative,  Sj;^(c,t)  are  small  at  t = ±1  and  know  also  that  S^(c,-t)  = 
S^(c,t)  for  even  order  and  S^(c,-t)  = -S^(c,t)  for  odd  order,  n.  One  can 
then  use  the  symmetry  property  to  find  S^(c,t)  in  the  interval  [-1,0] 
starting  with  a guess  for  S^(c,t)  at  t = -1,  say  S^(c,-1)  = 10"^.  Then 

s;(c,-l)  = S^(c,-1) 

Then  S^(c,t)  in  the  interval  [0,+l]  is  obtained  by  symmetry.  And  to 
obtain  S^(c,t)  in  the  interval  [+l,tg]  the  same  method  Runge-Kutta  and 
Adams-Moulton  is  used  with  starting  value  S^(c,+1)  and  S^;^(c,+1). 

6.4.2  Runge-Kutta  and  Explicit  Adams-Moul ton  Method 

The  second  order  differential  equation  for  the  prolate  spheroidal 
function  is  rewritten  as  in  Section  6.2.1  in  the  form  of  a system  of  two 
first  order  differential  equations. 


du, 

= f^  (t,u-]  *^2)  = U2 


du 


_2 

dt 


= f2(t,u^,U2)  = 


1 - t' 


'1  ^ 


(6.35) 


with  u-|  (-1 ) = guess  . . 

U2(-l)  = ( ^ 2 ^ ) • guess  . 

The  following  Runge-Kutta  (order  four)  is  used  to  find  the  starting  points 
for  Adams-Moulton. 

for  j = 0,1 ,2, ... ,N 


. -1  + jh 


where  h = (1/N)  [17] 
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and  = h 

k]2  = h 

^21  = ^ + (k^^/2),U2j  + (k^2/2)) 

k22  = h f2(t^.  + (h/2),u^j  + (k^T/2),U2j  + 

l<3i  = h f^(tj  + (h/2),u^j.  + (k2i/2),U2j  + (k22/2)) 

k32  = h f2(tj  + (h/2),u^j  + (l<2-|/2) ,U2j  + (K22/2)) 

*^41  " ^ '^31’^2j  42^ 

^42  ~ ^ "*"  ^31  ’ ^2j  "*"  ^32^ 

The  solution  at  the  next  mesh  point  j is  given  by 

'^l,j+l  " ‘^l.j  ^ ^*^21  ^4l  '^41^ 

and  U2  j^i  ~ ^2  j ^^12  ^^22  ^^32  "*"  ^42^ 


Then  an  explicit  Adams-fioulton  method  is  used  for  the  remaining  mesh 
points  as  follows 


^l,j+l  " ‘^l.j  ^^^2,j+l  ■ ^^2,j-l  ^2,j-2^ 


cHj^l  - X 


2t,. 


^2,j+l  ^2,j  ’■^^1  1 -^^+1  “2,j+l^ 


c^t?  - X 2t. 

1 - tj  ^Ij  1-  t?  ^2j) 
cH?  , - X 2t.  , 

* 1 '^j-1  1 - t]_i  ^2,j-l^ 

c^t?  2 - X 2t.  2 

^ ^ 1 - tj_2  ^l’j-1  1 - t?_2  ^2,j-2^  ^ 
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Rewriting  the  equations. 


“ij+l  ' “i.O  * * '^“2,0  ■ *“2,0-1  * “2,0-2^ 

“2,j+l  ' “2,0  **'^^''*  I-**®  “l,0+l  * * “2,0+1 ' 

+19(  eu,_j  +dU2_j) 

- “l,j-l  “2,0-1 > 


+ (g  Ui,j.2+  k l2,j-2)] 


a = 


0 = 


cH]^l  - X 
cH|  - X 


cH? 


ill 


^ - ^J-1 


cH,? 


-2  - X 


b = 


d = 


f = 


h = 


'^J±L 


2t. 

'1  _'^t'^ — 

^^j-2 
' - ^]-2 


bringing  the  expression  for  u-j  into  that  of  U2  j+-|  yields: 

“2,0+1  ° “2,0  * (•’''24)  [9(a{u,_j  +^  **“2,0+1  * '*“2,0  ■ *“2,0-1 

+ + b + 19(eu,_j  + d U2_j) 


- 5(e  + f + (g  u,_j_2+  k Uj.j.j)] 
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Separating  terms  and  solving  for  U2  j^.-|  gives 


(1  + 9 X 19  X a X (h/24)^  - 19  x d x(h/24))  u-  . 

^ 5 J 


+ (-9  X 5 X a X (h/24)2  - 5 x f x (h/24))  u,  . , 

^ j J “ I 


+ (9  X a X (h/24)^  + k x (h/24))  u 


2,j-2 


+ (9  X a X (h/24)  + 19  x e x (h/24))  u,  . 

' 9%. 


+ (-5  X e X (h/24)) 

+ (g  X (h/24))  j_2  ~\ 


(1  - 9""  X (h/24)2  X a - 9 X b X (h/24)) 


Now  that  u«  . , can  be  computed  explicitly  in  terms  of  u,  u«  . ,, 

^2  j-2  ^1  j’  *^1  j-1  ^1  j-2’  possible  to  find  an  explicit 

solution  for  u^  by  replacing  U2  j+-]  by  its  expression  to  get: 

“l,j+l  " ^l.j  a(h/24)"  - 19d(h/24))  U2^j 

- (45  a(h/24)^  + 5f(h/24))  U2^j_-, 

+ (9  a(h/24)2  + k(h/24))  U2^j_2 


l,j 


+ (9  a(h/24)  + 19e(h/24))  u 
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+ (-5  e(h/24)  u 


l.j-1 


+ g(h/24) 


(1  - 81(h/24)2  a - 9b(h/24))  + 19(h/24)  u 


2,j 


- U2_j.,  ^ (h/24)  U2_^.2 


This  is  the  explicit  Adams -Moulton  solution  for  u,  in  terms  of  (u-, 

I > j'*’ I I » j 

^l.j-1’  ^l,j-2’  ^2,j’  ‘^2,j-r  ^2.j-2)‘ 

Clearly,  the  solution  is  obtained  using  a Runge-Kutta  (order  four) 

with  a guess  for  u-|  j_2  and  U2  generate  the  solution  at  points  (j-1) 

and  (j)  and  then  use  the  three  step  explicit  Adams-Moulton  solution  for 

points  (j+l),(j+2)  and  so  on. 

Expansion  of  f(t)  in  the  interval  (-1.+1).  In  a first  step,  the 
eigenvalues  x^(c)  were  found  using  a numerical  method  (Runge-Kutta  and 
explicit  Adams-Moulton  method).  This  is  performed  by  program  EM.FTN  and 
the  result  is  in  the  table  on  the  following  pages.  Then  the  angular 
prolate  spheroidal  functions  in  the  interval  (0,  + t^)  were  found  using 
the  same  method.  And  the  expansion  is  performed  as  follows  for  a given 
bandlimited  function  f(t). 


Let  f(t)  = I Sp  <i>n(c,t) 


(for  all  t) 


to  simplify  the  notation  the  space-bandwidth  parameter,  c,  is  omitted. 
With 

♦-(t)  ' IT- 


n 


Table  1 
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Eigenvalues  for  c = 50  and  n = 0,1,..., 75. 


n 

Eigenvalue 

0 

.49246000 

2 

1 

.14823000 

3 

2 

.24619000 

3 

3 

.34311000 

3 

4 

.43897240 

3 

5 

.53375720 

3 

6 

.62744480 

3 

7 

.72001520 

3 

8 

.81144480 

3 

9 

.90171120 

3 

10 

.99078720 

3 

11 

.10786460 

4 

12 

.11652580 

4 

13 

.12505890 

4 

14 

.13346060 

4 

15 

.14172690 

4 

16 

.14985370 

4 

17 

.15783610 

4 

18 

.16566910 

4 

19 

.17334680 

4 

20 

.18086240 

4 

21 

.18820820 

4 

22 

.19537530 

4 

23 

.20235320 

4 

24 

.20912890 

4 

25 

.21568640 

4 

26 

.22200470 

4 

27 

.22805400 

4 

28 

.23378960 

4 

29 

.23914310 

4 

30 

.24404040 

4 

31 

.24851900 

4 

32 

.25287310 

4 

33 

.25751120 

4 

34 

.26264200 

4 

35 

.26825400 

4 

36 

.27427200 

4 

37 

.28064150 

4 

38 

.28732770 

4 

39 

.29430800 

4 

40 

.30156650 

4 
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n 

Eiqenvalue 

41 

.30919190 

4 

42 

.31687430 

4 

43 

.32490660 

4 

44 

.33318320 

4 

45 

.34169780 

4 

46 

.35044680 

4 

47 

.35942530 

4 

48 

.36862990 

4 

49 

.37805710 

4 

50 

.38770350 

4 

51 

.39756520 

4 

52 

.40763930 

4 

53 

.41792260 

4 

54 

.42841110 

4 

55 

.43910230 

4 

56 

.44999350 

4 

57 

.46108150 

4 

58 

.47236470 

4 

59 

.48384180 

4 

60 

.49551320 

4 

61 

.50738030 

4 

62 

.51944430 

4 

63 

.53171110 

4 

64 

.54418530 

4 

65 

.55687320 

4 

66 

.56978160 

4 

67 

.58291700 

4 

68 

.59629310 

4 

69 

.60991060 

4 

70 

.62378010 

4 

71 

.63790870 

4 

72 

.65229760 

4 

73 

.66695260 

4 

74 

.68188000 

4 

75 

.69708020 

4 
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where 


= 


+1 


-1 


S^(t)  dt 


and 


n " A. 


+1 


f^it)  cj)^(t)  dt 


+1 


-1 


it 

n 


a_  = 


1 


" u A" 
n n 


+1 


-1 


• s^(t)  dt 


where  fg(t)  is  the  observed  data  (or  measured).  The  expanded  function 
f(t)  win  be  given  by 


+1 


J -1 


f(t)  • S (t)  dt 


f(t)  = I 

n u_  A. 


f(t) 


+1 


-1 


n n 


f(t)  S^(t)  dt 


[-T  S„(t)l  (6-36) 


n u' 


s„(t)  . 


(6.37) 


This  is  the  expression  of  f(t)  in  the  interval  -1  < t < +1  in  a series  of 
prolate  spheroidal  functions. 

Extrapolation  in  the  interval  [+1,°°].  In  order  to  evaluate  Equation 
6.37  in  the  region  t>l,  properly  normalized  S^(t)  must  be  available  as 
well  as  values  for  u^.  Here  we  will  obtain  the  normalized  S^(t)  by  con- 
tinuing the  numerically  calculated  S^(t)  with  an  approximate  analytic 
form  valid  for  large  t (t  > t^).  Rewriting  the  integral 
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+ 00 


S^(t)  dt 


= 2 


S^(t)  dt 


= 2 [ 


S^(t)  dt  + 


S^(t)  dt  + 


S^(t)  dt  ] 


Let 


uf  = 2 


S^(t)  dt 


v!  = 2 


S^(t)  dt 


w!  = 2 


S^(t)  dt 


and 


n n n n 


Now,  for  t^  large,  one  approximates  the  angular  prolate  spheroidal  func- 
tion S^(t)  by  the  radial  function  [12]  in  the  form  of  a cosine  function 


^ cos  (ct  - (n  + 1)  tt) 


so  that 


wj  = 2 


cos^Ect  - -^{n  + 1)tt]  dt 


set 


a = c 


b = --j  (n  + 1 )ir  X = t 


to  write 


95 


cos  ^(ct 


(ax): 


and 


cos^(ax  + 


[ Y (1  ■*■  cos  2(ax  + b))] 


cos^(ax  + b)  _ 1 ^ cos  2b  cos  2ax  sin  2b  sin  2ax 

(ax)  2 ■ 2(ax)^  2a ^ ' x^  2P"  ' “7^ 


and 


cos^(ax  + b)  j..  _ 1 


1 ^ cos  2b  r cos  2ax  o 
7 ~2i^  1- 2a 


sin  2ax 


dx] 


sin  2b  r sin  2ax  , o. 
— 2p-  [ + 2a 


cos  2ax 


dx  ] 


with  b = - (n+1)iT  n = 0,1,2,... 


^0  " ■ 2 


b.j  = -7T 


b2  - 2 ^ 


2b  = -IT 


■2n 


-3v 


Thus  sin  2b  = 0 


And 


cos^(ax  + b) 


dx  = 


1 , cos  2b  f cos 

7a^  ^ ~2P“  T 


2b 


- 2a 


sin  2ax 


dx) 


tr 


2^  2a^  - a 


sin  2ax 


dx 


0 Jo 
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r" 


Let  I 


1 


sin  2 ax  j..  1 

^ ' 2 


+00 


sin  2ax 


dx 


I]  = Y [I]  where  I = 


„i2az 
e 


dz 


Using  a contour  integral,  I is  found 


and 


I = iir 


I = - 
M 2 


Therefore 


cos^jax  + b)  , ^ 1 (-1)”'*'^  ^0 


(ax): 


2a^  2^ 


n+1 


and 


"n  = 2 


cos^(ax  + b) 
ax)^ 


dx 


1 / , Nn  cos  2a  t 


0 + 


IT 

I 


sin  2ax 


dx  ] 


.-2 

a 


sin  2ax 


dx 


where 


I = 


sin  2ax 


dx  = 


2at, 


dz  = Si  (2atQ) 


can  be  evaluated  numerically  using  a Simpson  rule,  and. is  tabulated  [18], 

T / 1 cos  2a  t„  / i sn 

21] 


n a^t. 


= u^  + 

n n n n 


and 


is  determined 
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with  for  the  interval  0 < t < +1 

vf  for  the  interval  +1  < t < t 
n 0 


for  the  interval  t„  < t < “ 
n 0 


Therefore,  outside  the  interval  -1  < t < +1,  the  extrapolation  of 
the  bandwidth  limited  function  f(t)  will  be  given  by  an  expression  similar 
to  Equation  6.37  where  the  coefficient  u^  in  the  denominator  is  replaced 


byA^. 


f(t)  = I 
n 


+1 


-1 


f(t)  S„(t)  dt 


n 


Sn(t) 


(6.38) 


This  scheme  was  applied  using  programs  NZA.FTN  and  NSDBLE.FTN. 

Another  method  of  continuing  the  S^(t)  to  “was  also  attempted:  After 
the  S^(t)  were  calculated  in  the  interval  +1  < t < +5,  they  were  extended 
using  the  integral  relationship 


s„(t) 


■+1 
• -1 


sin  c(t  - s) 
TT(t  - s) 


Spts) 


ds 


[12] 


At  this  point,  the  angular  prolate  spheroidal  functions  for  the  value 
of  the  space- bandwidth  parameter  c = 50  were  now  known  in  the  interval 
-5  < t < +5  and  they  were  saved  on  disk. 

Beyond  t = 5 these  functions  are  approximated  by  the  expression 


S^(t)  = AM 


n 


cos(ct  + 6) 
ct 


(6.39) 


and  its  first  derivative 
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S'(t)  = AM  [~c  sin(ct  + 6)  • ct  - c cos(ct  + 9)  j 


S"(t)  = - AM  [ sin(ct  + 6)  ^ cos(ct  + 6) 

" " t ct^ 


(6.40) 


The  angular  prolate  spheroidal  function  and  its  derivative  were 
evaluated  at  t = 5 

Sn(5)  = AM^  (6.41) 

and 

S'(5)  = - AM„  [ ^ ^ ^ ] (6.42) 

The  ratio  of  the  two  functions 
Sn(5) 

RS  = = f(sin  9,  cos  9) 

Sp(5)  was  known  and  S^(5  - At)  was  also  known.  And  using  a backward 
difference  operator  S^(5)  is  known,  too.  Therefore  the  ratio,  RS,  is 
known.  Hence,  using  a numerical  approach  one  can  determine  the  value 
of  the  parameter  9 that  yields  the  same  ratio,  RS.  And  replacing  9 by 
its  value  in  Equation  6.41  will  give  the  amplitude,  AM^,  of  the  angular 
function  for  t greater  than  5,  and  consequently,  S^(t)  will  be  known  for 
any  t. 

The  next  step  was  to  normalize  these  functions  using 

f 

S^(t)  dt  = 1 

J — oo 

This  will  determine  the  value  of  the  normalization  constant  A to  obtain 

n 

the  prolate  spheroidal  functions  4>p(t)  = A^S^(t)  and  a bandlimited  function 
f(t)  as  follows: 
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= I bp  4>p(t) 

n 

where 

C+1 

f(t)  (J)^(t)  dt 

•-1 


4>p"(t)  dt 

•-1 


(6.43) 


(6.44) 


6.5  Conclusion 

The  above  scheme  was  applied  to  extrapolate  two  functions.  One  was 
a cosine  wave  in  -1  < t < +1,  and  the  other  was  a sinc-function.  The 
expansion  in  the  interval  -1  < t < +1  was  perfect  but  the  extrapolation 
was  not  as  expected;  however,  it  looks  better  for  the  sinc-function  case 
than  for  the  cosine- function  case  as  depicted  in  Figure  6.3  and  Figure  6.3, 
respectively. 

The  programs  used  were  ANTHET.FTN,  CORFAC.FTN,  EXXDl.FTN,  and  EXXD2.FTN. 

We  believe  that  the  continuation  processes  used  for  t < 5 introduced 
significant  errors  in  the  normalization  for  large  order  (large  n)  eigen- 
functions. The  resulting  eigenvalues  which  range  from  Xq=1  to  X^q^IO"^®, 
for  example,  caused  the  extrapolations  to  fail.  As  we  will  see  in  the 
next  chapter,  more  accurate  calculations  seem  possible  in  the  frequency 
domain  rather  than  the  time  domain. 


TOO 
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CHAPTER  VII 

PROLATE  SPHEROIDAL  METHODS  IN  THE  FREQUENCY  DOMAIN 


The  methods  proposed  in  the  chapter  have  not  yet  been  tested  and  are 
included  as  suggestions  for  future  work.  They  are  an  outgrowth  of  the 
work  described  in  Chapter  VI  which  showed  the  difficulties  involved  in 
obtaining  and  manipulating  prolate  spheroidal  functions  with  large  values 
of  c (time  limit-band  limit  product). 

We  describe  below  a method  which  may  be  used  to  obtain  a function 
g(w)  which  has  the  properties: 

i)  when  convoluted  with  the  observed  Fourier  spectrum,  it  narrows 
the  resolution  function. 

ii)  It  does  not  require  the  use  of  very  high  order  prolate  spheroidal 
functions 

iii)  The  resulting  resolution  function  is  not  necessarily  a sine  (u) 
but  may  be  chosen  by  the  analyst  in  order  to  reduce  the  effects 
of  random  noise. 

In  essence  our  original  problem  was  to  find  a function  g(o))  such  that 
when  operating  on  a sine  function  from  to  +«>  would  yield  a 6-function  or 
a narrower  sine  function.  It  was  believed  that  this  operation  might  be  a 
convol ution: 

,+co 

g((jj)  ■ sine  t(u)q  - u)  do)  = (7.1) 

J — OO 

or  = narrower  sine  function 
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or,  for  finite  data  sets,  is  it  possible  to  find  g(  ) such  that 

g((i))  : sine  t(cjo^  - co)  doj  = 6(w^)  (7.2) 

• -n 

or  = narrow  sine  function? 


Replacing  in  Equation  7.1  and  Equation  7.2  the  sine  t(o)^  - u)  and 
the  by  their  expression  in  terms  of  the  prolate  spheroidal  functions 

[11]  given  by 


sine  - 0))  = ^ ^ 


(7.3) 


“ KM  4>_(o)) 
6(„  - “)  = I J 

On 


(7.4) 


Equation  7.2  becomes 


■K2 


“ “ (b  (o)  ) (b  (o) 

gM  [ ^ I tJu.)]  doj  = I ” °.  — 

T 0 n 0 n Q 


(7.5) 


of  course 


g(cj)  = 0 |wl  > (bandlimited) 


Rearranging  terms  in  Equation  7.5, 

r+fl 


- 1 4'  (w  ) [ 

T g ^n^  o' 


~n 


<t>n(w  ) 4>_(o) 

g((o)  (^„(o))  do)  ] = 

0 n 


(7.6) 


Let  g(ijj) 
with 


= 


+00 


g(oj)  4>^(oj)  dw 


4^2 


g(‘^) 


4f,(w)  dto 


(7.7) 


Using  Equation  7.7  and  replacing  in  Equation  7.6  the  expression  between 
the  brackets  yields: 


00 


00 


giving  for  the  coefficient  a 


n 


a 


n 


And  the  sought  after  function  will  be 


(7.8) 


Clearly  the  determination  of  the  above  derived  function  is  subject 
to  a truncation  of  the  summation  computed  up  to  an  order  N instead  of 
infinity.  We  find  that  the  function  g(oj)  is  a sine  function  for  values 
of  N up  to  a critical  value  (2c/7t),  and  above  this,  due  to  the  term 
in  the  denominator  of  Equation  7.8,  the  higher  orders  predominate,  the 
sine  function  is  destroyed,  and  the  function  g(o))  takes  on  the  shape  of 
the  higher  order  prolate  spheroidal  functions,  What  is  happening 

seems  to  be  that  the  functions  g(o))  requires  high  order  (hign  n)  in  order 
to  describe  it  accurately  and  we  were  able  to  compute  these  high  order 
functions  sufficiently  accurately. 

However,  when  using  Equation  7.4  to  calculate  a 6-function,  we  obtain 
a sine  function  20  points  wide  for  a value  of  N = 30  and  a narrower  sine 
function  10  points  wide  for  N = 75.  As  one  gets  much  closer  to  a 5-function, 
say,  one  or  two  points  wide  for  a very  large  number  of  orders,  this  resembles 
the  deconvolution  function,  g(o)),  one  point  wide  obtained  by  direct  inver- 
sion of  the  convolution  matrix,  R,  (Chapter  III)  constructed  with  a sine 
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function  as  resolution  function, 

E.  ■ £ “ A » 

as  depicted  in  Figure  7.1. 

Therefore,  two  points  must  be  emphasized,  the  first  one,  no  function 

g(co)  can  be  found  which,  when  convoluted,  i.e., 

+ 00 

g(o))  sine  T(a)^  - to)  dco 

• . 00 

with  a sine  function  produces  a 5-function. 

The  second  point,  if  the  convolution  is  truncated,  i.e., 

g(to)  sine  T(to^  - to)  dto  , 

then  such  a g(to)  can  be  found  as  in  Equation  7.8,  but  it  requires  high 
order  eigenfunctions  for  accurate  computation. 

If  one  is  content  with  partial  deconvolution,  i.e.,  turning  the 
sine  t(o))  resolution  function  into  a narrower  resolution  function  (though 
not  a delta  function),  such  high  order  eigenfunctions  may  not  be  needed. 

(It  is  not  necessary  to  extrapolate  f(t)  to  extrapolation  from  t^  to 
5t^,  for  example,  will  narrow  the  sine  (xw)  by  a factor  of  5). 

This  reveals  one  additional  exciting  possibility:  Sine  the  new 

resolution  function  shape  may  be  chosen  by  the  analyst,  it  may  be  chosen 
in  order  to  reduce  the  effects  of  random  noise.  Let  us  examine  this  possi- 
bility. The  observed  spectrum  is  the  convolution  of  the  true  spectrum  with 
a sine  function. 


''obs<“>  ' ''true  “ 


We  know  that  one  can  find  a deconvolution  function  g(o))  such  that 
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C 

o 


t/> 

c 

o 


s- 

3 
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g(uj)  B sine  (tco)  = 6((d) 

as  expressed  in  Equation  7.2.  But  this  is  very  difficult  to  calculate 
because  we  need  high  order  eigenfunctions.  So  for  less  than  perfect 
deconvolution  we  want  to  find  g^(o3)  such  that 

g"(o))  a sine  (tw)  = Ga(co) 


where  Ga(w)  is  the  new  resolution  function  chosen  by  the  analyst  as  a 
Hanning,  a Hamming,  a Bartlett,  or  a Gaussian  window  for  example. 

Let  us  suppose  that  Ga(co)  is  a Gaussian 

N 

Ga(oj)  = J a (j)  (w) 
n=0  " " 

where  only  the  first  30  orders  are  needed,  N = 30,  and  let  us  find  g"(u3). 
Let 

CO 

g"(u)  = I b (J)  (oi) 

0 n n 


and  we  wish  to  satisfy 


r+Jl 


g"(w)  sine  x(y  - cj)  = Ga(y)  = I 


Ga(ii))  du) 


r+Jl 


g"(o))  4>^(u)  do) 


a^  0 n > N 

Multiplying  through  by  4>j^(y)  and  integrating, 

•+f2 

g"(w)  ^ (j)^(a))  du)  = a^ 

■ 
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and 


00 


= I 
0 


and  sine  a_  = 0 


n > N, 


n 


g''(w)  = I ^ <}>„(w)  and  only  the  first 

0 ^ ^n  " 


30  orders  are  needed  to  achieve  a partial  deconvolution  and  better  resolution. 

One  may  say,  roughly,  that  if  Ga(u)  is  five  times  narrower  than 
sine  (xto),  the  above  calculation  is  equivalent  to 

(1)  extrapolating  f(t)  by  a factor  of  5, 

(2)  windowing  f(t)  by  the  inverse  Fourier  transform  of  Ga(o)),  and 

(3)  Fourier  transforming  the  result  to  obtain  an  improved  Fourier 


spectrum.  It  is  important  to  note  that  the  windowing  function 
may  be  chosen  by  the  analyst  in  order  to  reduce  the  effects  of 
random  noise  on  the  improved  spectrum. 


CHAPTER  VIII 
CONCLUSIONS 

Four  methods  of  deconvolution  of  a Fourier  spectrum  have  been  studied 
here:  the  inverse  matrix  method,  the  nibbling  technique,  the  iterative 

technique,  and  the  prolate  spheroidal  functions  method.  General  conclusions 
are  as  follows: 

(1)  It  has  been  observed  that  the  matrix  methods  for  deconvolution 
were  of  limited  value,  and  became  extremely  sensitive  tho  noise 
when  the  matrix  size  increased  to  correspond  to  real  life  problems. 
A recursive  method  of  solution  was  described,  however,  which 
allows  rapid  inversion  of  very  large  matrices  of  Toeplitz  type. 

(2)  The  nibbling  technique  is  quite  promising,  particularly  in  the 
case  of  widely  spaced  peaks  in  the  Fourier  spectrum.  It  was  ob- 
served that  when  the  peaks  came  closer  and  closer  to  each  other 
(less  than  one  main  lobe  apart)  there  was  creation  of  satellite 

or  spurious  peaks  around  the  main  peaks.  However,  instead  of  just 
nibbling  away  the  observed  Fourier  spectrum  with  a resolution 
function  located  at  the  occurrence  of  the  maximum  (equivalent  to 
cross-correlating  the  observed  spectrum  with  a one-point  response 
function),  one  could  use  the  following  procedure: 

(i)  Cross-correlate  the  observed  spectrum  F(w),  with  an  M-point 
resolution  function,  R(w)  (perhaps  just  the  N points  which 
form  the  main  lobe), 

(ii)  Locate  the  maximum  of  the  cross-correlation  function. 
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no 


(iii)  Take  away  a nibble  from  the  Fourier  spectrum  and  repeat 
the  procedure  until  either  of  e,  a,  or  M (defined  in 
Chapter  IV)  is  reached.  For  closely  located  peaks  this 
procedure  is  promising  and  will  avoid  the  creation  of 
satellite  peaks  in  the  deconvoluted  spectrum.  In  this 
manner  the  spectral  resolution  will  be  improved,  but  at 
the  expense  of  a much  longer  computation  time. 

(3)  In  the  case  of  the  iterative  technique,  our  investigation  showed 
that  the  method  did  not  always  converge  in  amplitude  to  the 
desired  function  even  in  the  absence  of  additive  noise.  Further- 
more, since  improvement  of  resolution  in  the  frequency  domain  is 
equivalent  to  extrapolation  in  the  time  domain,  the  presence  of 
such  noise  would  greatly  affect  the  uncertainties  in  the  time 
extrapolation. 

(4)  The  prolate  spheroidal  functions  method  proved  excellent  for  the 
recovery  of  the  function  in  the  time  domain  as  well  as  in  the  fre- 
quency domain;  but  the  extrapolation  was  not  as  good  in  either  one 
of  the  domains  because  higher  order  eigenfunctions  were  needed  and 
it  was  very  difficult  to  obtain  higher  order  eigenvalues  with  suf- 
ficient accuracy.  A very  promising  alternative  was  noticed  however. 
Instead  of  time  extrapolation  as  first  proposed  by  Slepian  and 
Pollack  [12],  frequency  interpolation  offers  greater  computational 
ease  as  well  as  the  likelihood  of  significant  reduction  in  the 
effects  of  random  noise.  The  procedure  is  as  follows: 

(i)  Obtain  a set  of  prolate  spheroidal  eigenfunctions 

in  the  frequency  domain  for  a value  of  c much  larger  than 
usually  needed  (e.g.,  c = 1000,  t = 1,  n=  1000  where  the 
bandlimit  on  f(t)  is  « 1000. 


(ii)  Then  the  function 


g((^) 


00 


I 

n=0 


1000 

IT  X 

n 


when  convoluted  with  ^^^^(03)  would  deconvolute 
into  would  be,  however,  computationally 

impossible  since  high  order  eigenfunctions  would  be 
needed.  (The  presence  of  random  noise  would  also  greatly 
increase  the  difficulties  since  this  would  be  equivalent 
to  extrapolating  f(t)  to  infinity— including  the  noise), 
(iii)  Instead,  one  may  decide  partially  to  deconvolute  such  that 


''obs^'^"^  " 8 sine  (xu)) 


true 


becomes 


^improved^^  ^ ~ ^true^'^^  ® Ga(o)) 


where  6a(u)  is  several  times  narrower  than  sine  (xw),  ■ 
but  not  so  narrow  that  high  order  eigenfunctions  are 
required  to  express  it.  Then  when 

low  N 

Ga(co"  - 0))  = I cf)_(a)) 

0 n n 

the  partial  deconvolution  function  g^u)  becomes 
0 " ""  ^n 

and  only  a small  number  of  terms  may  be  necessary,  and 
the  and  41^(0))  will  be  tractable. 

A significant  advantage  to  this  method  of  interpolation  in  the 
frequency  domain  is  that  the  final  resolution  function  is  chosen 
by  the  analyst  (it  need  not  be  a sine  function)  and  it  may  be 
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chosen  to  minimize  random  noise  effects.  It  is  equivalent  to 
extrapolating  in  the  time  domain  while  simultaneously  windowing 
or  tapering  the  extrapolation,  thus  reducing  the  effects  of 
extrapolating  random  noise. 


APPENDIX  A 

ORGANIGRAMS  FOR  DECONVOLUTION  WITH  MATRIX  METHOD 


A.1  Deconvolution  with  Matrix -I 


B = 
BB  = 

Me  = 
N = 


> ( 


row  of  A"^ 

working  vector 

observed  spectrum 

estimated  spectrum 

size  of  the 
deconvolution  matrix 


1 

C 
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1]5 


B C 


Figure  A.l  Deconvolution  with  inverse  matrix  read  one  row  at  a time 
with  REWIND  operation. 
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A. 2 Deconvolution  with  Matrix -I I 
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Figure  A. 2 Organigram  to  read  the  deconvolution  matrix  one  row  at  a 
time  without  REWIND  operation. 
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A. 3 Implementation  of  the  Algorithm  in  Figure  A 


1=1 


KZ(1)  = 

1 

J = 

= 1 

■ 

KZ(J)  = 

1 

LOC(l)  = 

1 

1 

2 

4 

7 

11 

1 

L0C(2)  = 

2 

2 

2 

3 

5 

8 

12 

3 

L0C(3)  = 

4 

4 

4 

5 

6 

9 

13 

5 

L0C(4)  = 

7 

6 

7 

8 

9 

ig 

14 

7 

L0C(5)  = 

11 

8 

11 

12 

13 

14 

Jl 

9 

10 

L0C(1)-1 

= 0 

11 

12 

read  nothing 

13 

14 

1 read  B(1 ) 

15 

J = 2 


J = 3 


J = 4 


J = 5 


(L0C(1)+1)  = 2;(L0C(2)-1)  = 1 
read  B(2) 

(L0C(2)+1)  = 3;(L0C(3)-1)  = 3 
read  'JUNK'  (element  3) 
read  B(3) 

(L0C(3)+1)  = 5;  (L0C(4)-1)  = 6 
read  'JUNK'  (elements  5 and  6) 
read  B(4) 

L0C(4)+1  = 8;(L0C(5)-1)  = 6 
read  'JUNK'  (elements  8,  9,  and  10) 
read  B(5) 


■2  for  the  Case  N=5 
3 4 5 

6 10  15 

elements  of 
row  1 are  read 
in  one  sweep 
without  rewinding. 
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1=2 


LOC(l)  = 2 
LOC(2)  = 3 
L0C(3)  = 5 


read  in  one 
sweep  of  the 
data. 


elements  of 
row  2 are 


L0C(4)  = 8 
L0C(5)  = 12 

LOC(l)  -1=1 

J = 1 read  'JUNK'  (element  1) 
read  B(1 ) 

J = 2 read  B(2) 

J = 3 L0C(2)  +1=4;  L0C(3)  -1=4 

read  'JUNK'  (element  4) 
read  B(3) 

J = 4 L0C(3)  +1=6;  L0C(4)  -1=7 

read  'JUNK'  (elements  6 and  7) 
read  B(4) 

J = 5 LOG  (4)  + 1 = 9 ; L0C(5)  - 1 = 11 

read  'JUNK'  (elements  9,  10  and  11) 
read  B(5) 


1 
11 
12 
TT 

14 

15 


c\j  |oo  1"^  IT)  |i£)  cx)  |cn  o 
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1=3 


LOC(l)  = 4 
L0C(2)  = 5 
L0C{3)  = 6 
L0C(4)  = 9 
L0C(5)  = 13 


entries  of 
row  3 are 
read  in  one 
sweep 


LOC(l)  -1=3 

read  'JUNK'  (elements  1,2,  and  3) 


11 

12 

13 

14 

15 


read  B(l),  B(2),  B(3) 

J = 4 L0C(3)  +1=7  ; L0C(4)-1  = 8 

read  'JUNK'  (elements  7 and  8) 
read  B(4) 

J = 5 L0C(4)  + 1 = 10  ; L0C(5)-1  = 12 

read  'JUNK'  (elements  10,  11,  and  12) 
read  B(5) 


cvj  CO  ^|ir)|uD|r^  oo  ch|o 


12J 


LOC(l)  = 7 

elements  of 

L0C(2)  = 8 

row  4 are 
read  in  one 

L0C(3)  = 9 

sweep 

L0C(4)  = 10 

L0C(5)  = 14 

LOC(l)  -1=6 

read  'JUNK'  (elements  1 through  6) 
read  B(l) 

B(2) 

B(3) 

B(4) 

J = 5 L0C(4)  + 1 = 11  ; L0C(5)  - 1 = 13 

read  'JUNK'  (elements  T1 , 12,  and  13) 
read  B(5) 
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LOC(l)  = 11 
LOC(2)  = 12 
L0C(3)  = 13 
L0C(4)  = 14 
L0C{5)  = 15 

LOC(l)  - 1 = 10 

read  'JUNK'  (elements  1 through  10) 
read 

B(l) 

B(2) 

B(3) 

B(4) 

B(5) 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

J2 

13 

14 

II 


A. 4 Storage  of  a Persymmetric  Matrix 


12.4 


Figure  A. 3 Organigram  for  the  location  of  an  element  of  a persymmetric 
matrix  stored  in  a vector  form. 
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A. 5 Implementation  of  the  Algorithm  in  Figure  A. 4 for  the  Case  N=6 


A Fortran  program  for  the  algorithm  of  Figure  A. 3 is  NEWLOC.FTN 
and  is  verified  in  the  following  example  of  a 6 x 6 matrix. 


N=6 


I J 


1 

1 

1 

1 

1 

1 


1 

2 

3 

4 

5 

6 


2 

2 

2 

2 

2 

2 

3 

3 

3 

3 

3 

3 


1 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 


1 

2 

3 

4 

S 

6 

2 

7 

8 

9 

20 

5 

3 

8 

11 

12 

9 

4 

4 

9 

12 

11 

8 

3 

5 

10 

9 

8 

7 

2 

6 

5 

4 

3 

2 

1 

K(I,J) 

1 

2 

3 

4 

5 

6 

2 

7 

8 
9 

10 

5 

3 
8 

11 

12 

9 

4 


All  the  underlined  elements 
(I,J)  correspond  to  their 
vector  location  K. 


A. 6 Deconvolution  with  Matrix-III 
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I+l 


Figure  A 


Read  A'^ 
Read  M 

0 

I 

I = 1,  N/2 

^ 1 


II  = N+l-I 

1 

J=l, 

N 

»■■■  I 


JJ  = N+l-J  ; K 

DM  = DM  + A"^(K).Mq(J) 
DN  = DN  + A'^(K)-Mq(JJ) 


Estimated 

spectrum 


4 Restoration  of  the  estimated  spectrum 
matrix  stored  in  a vector  of  N(N+2)/4 
storage). 


= deconvolution  matrix 
= observed  spectrum 
= estimated  spectrum 


with  deconvolution 
locations  (persymmetric 
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A. 7 Recursive  Solution  of  Simultaneous  Equations 
Involving  a Persymmetric  Matrix 

A. 7.1  Presentation 

Consider  the  system  of  equations 


where  is  a square  persymmetric  matrix  with  entries  rg,  r-j , >"2”‘"’’"m 
assumed  known,  and  is  a vector  with  components  gg,  g-j , g2,...,gj^  which 
is  also  assumed  known.  At  a given  step  n of  the  Toeplitz  recursion,  the 
prediction  error  coefficients  satisfy  the  equation  [ig] 


^®n0’  ^nl’- •*’^n,n-r  ^nn^  ^ ^“n’ 

0,  0,  0,. . . ,0) 

(A.2) 

and 

^^nO’  ^nl '^n,n-r  ’*'nn^  ^^0’ 

9l , • • • , 9p ) 

(A.3) 

At  the  completion  of  step  n,  the  following  terms  are  stored 

^nO’  ^nl ®nn  ’ “n’  ^n 

and 

^nO’  ■'"nl ^nn  ’ ^n 

such 

that 

^^nO’^nl ^nn’  Vl  " ^“n’ 

0,  0,...,0,  3^) 

(A.4) 

and 

^^nO*'''nl’'-'’^nn’  Vl  " 

9l’---’  ^n^ 

(A.5) 

R^^-j  is  a Toeplitz  persymmetric  matrix.  Equation  A. 4 can  be  written  as 
follows 


(0,  a 


nn’  '"n,n-l ‘‘nl’ 


®n0> 


'n+1 


0,  0, . , 


(A. 6) 


12? 


Multiplying  Equation  A. 6 by  k 


"n  "n  %,n-l =nT  ''n  "‘nO*  Vl  ' 


(k^  6^,  0,  0,...,  0,  k^  a^)  (A. 7) 


Adding  Equations  A. 4 and  A. 7 


^^nO’  ^nl  *^n  ^nn’"*’  ^nn  “^n  ®nl’  *^n  ®n0^  ^n+1 


(“n  ■*■  '^n  ^n’  ^n  '''  *^n  “n^ 


It  is  desired  that  Equation  A. 8 be  equivalent  to 


^®n+l,0’  ^n+l,r---’  ®n+l,n*  ®n+l,n+l^  Vl 


(cx^_j_1  j Oj  0,...t0) 

Comparing  Equations  A. 8 and  A. 9 yields 

^n+1,0  " ®n0 

n+1,1  nl  n nn 


a„.,  „ ' = a„„  + •<„  a„, 

n+l,n  nn  n ni 


and 


^n+l,n+l  ^n  ^nO 


+ k a = 0 »-  k„ 

n n n n 


a = a + k 3 
n+1  n n n 


6 

a 


n 

n 


(A.9) 


(A. 10) 


(A.ll) 


(A.12) 
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Equation  A.4  ► Vl.ntV  “>  W = 


» 0,  0,  •••>  0,  ^p+] ) (A.  13) 

Equation  A. 13  ^ + ...  + 

(A.14) 

Using  Equation  A. 13,  the  recursion  step  (n+1)  is  ready  to  be  started. 
Rewriting  Equation  A. 9 as  follows 

^®n+l,n+l’  ^n+l,n’“”^n+l,l’  ^n+1,0^  ^n+1  ^ O.---.  0,  a^^.j ) 

(A.15) 

Multiplying  Equation  A.15  by  q^, 

^^n  ^n+l,n+r  ^n  ^n+l.n’"-’  ^n  ^n+1,1’  ^n  ^n+1,0^  ^n+1  " 


(0,  0,...,  0,  q^  (A. 16) 

Adding  Equations  A. 5 and  A. 16, 

*^n0  * ‘'n  ®n+l,n+l’  *"01  *"nn  ‘'n  ®n+l,l’  % ^n+1,0* 

•''n+1  = <30’  9n-  % * % Vl> 


It  is  desired  that  Equation  A. 17  be  equivalent  to 

^Vl,0*  "^n+l,!’--- ’ Vl,n+1^  .jTjl  " ^9q’  5r-”’Vl^ 


Comparison  of  Equations  A. 17  and  A. 18  yields 
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and 

and 

Hence 

Assuming 

Equation 


Vi,o 

^n,0  ^n  ^n+1  ,n+1 

^n+l , 1 

= ^n  ^n+l,n 

• 

(A. 19) 

^n+1 ,n 

" ■''nn  ^n  ^n+1,1 

^n+1 ,n+l 

" ”n  ®n+l,0 

^n+l  " 

1 ^n  '^n+l  ’ *^n  ^^n+1 

^n>/Vl 

(A. 20) 

Equation 

A. 5 ^ Vi,0  V2 

Vi,i  Vi 

Vl,n+1  ’"l 

(A. 21) 

Vl,r--”Vl,n+l’  '^n+2 

= 

(gg*  9i>  gp>  gp+i*  Yp+i)  (a. 22) 


^ ^ “0  ^00  ’ ’"o  ^ “0  ~ *"0 

9g  is  known 

®oo  '"0 " “0 


Equation  A. 4 
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Equation  A.5 . f = g^/r^ 

^00  ’"o  " 9q 

^00  "'i  = ^0  ^0  = ( V'^q)  • "'l 

®00  " ^ 

»o  = >"0  ^'"o  S^ven) 

9o  ( given) 

6q  = r^  (r^  given) 

^00 " 9o^'"o 
^0  " "^00*  *"1 

A. 7. 2 Summary  of  the  Toeplitz  Recursion  Described  Above 
Given  the  convolution  matrix 

(tq,  r^,  r^) 

and  the  vector 

^2 Vr 

at  the  completion  of  step  n, 

®n0’  ®nr  ®n2”--’  ®n,n-l’  ®nn 
and 

^nO’  "^nr  ‘^n2’--*’  "^n.n-l’  "^nn 

are  known.  To  execute  step  (n+1), 

(1)  Find  with  Equation  A. 11, 

(2)  Find  ®n+l,n+l 


“n’ 


T, 


with  Equation  A. 10 
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(3)  Find  and  with  Equations  A.  12  and  A.  14, 

(4)  Find  with  Equation  A. 20. 

(5)  Find  the  set  '*'n+l,l ^n+l,n+l^  Equation  A. 19. 

(6)  Find  Yp+-|  with  Equaton  A. 21. 

(7)  Increment  (n+1)  to  (n+2)  and  go  to  Step  1. 


The  final  values  of  f ^ f o.-.-.  are  the  solutions  of  the 

mO  m I m2  ran 

system  of  Equation  A.l.  Starting  with  n = 0 and  the  values  of 

^00  ^ ^0 " '"o’  ^0 ""  ""r  ^00  ^ ^^o^'"o^  ' ""i 

recursively  one  ultimately  obtains  the  solution  of  the  system  of  equations. 

A. 7. 3 Example  of  Using  the  Toeplitz  Recursion  Algorithm 
Let 

Q “1  » ^1  “ .75  , r 2 ~ .5  , r^  ~ .25  » r^  ~ 0 

9g  ~ 0 , g-|  - , 25  , ^2  ~ ^ ’ ^3  ~ ’ 9^  ~ 0 

Step  n = 0 

®00  " ^ 

°^0  ^ '"o 

Bo  = r^  = .75 

^00  ^ ^o'^’"o  ^ ° 

Yo  = (gg/'^o)-''!  = ° 

Step  n+1  = 1 

*^0  ^ “ ^^0^“0^  = -.75 
^1,0  " ®00  " ^ 

a.|1  “ kQ  agg  — — .75  x 1 ~ *,75 
a.,  = ap  + kp  Bg  = 1 + (-.75)(.75)  = .4375 
Bi  = a.|Q  r2  + a.j.|  r.j  = 1 x (.5)  + (-.75)(.75)  = -.0625 
90  = (9i  - Yg)/a^  = (.25  - 0)/(.4375)  = .5714 
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f^O  = -^00  ^0  ° (.5714)- (-.75)  = -.4286 

f.,^  = qp  a^Q  = (.5714)(1)  = .5714 

Yi  = f.]g  i"2  + f-|i  r-|  = -.4286  x .5  + .5714  x .75  = .2143 
step  (n+1 ) = 2 

= -(6^/a^)  = - (. 0625)7.4375)  = .1429 

®20  ^ ®10  " ^ 

^21  " ®11  “^l  ®11  " (.1429)(-.75)  = -.8572 

®22  " *^1  ^10  " ^ 

= .4375  + .1429  (-.0625)  = .4286 

h " ®20  *"3  ^21  ’"2  ®22  ^1  " ''•(•25)  + (-.8572)(.5)  + (.1429)(.75) 

= -.0714 

°1  " (^2  - Y-,)/a2  = (1  - .2143)7.4286  = 1.8332 

f2o  = fio  ^ ^1  ®22  ^ --4286  + 1.8332  x .1429  = -.1666 
^21  " ■''ll  ^1  ®21  " 1.8332  X (-.8572)  = -1 

f22  = 320  ■'•8332  x 1 = 1.8332 

^2  ” ^20  ’"3  ^21  '^Z  ^22  *^1  “ -•1666  x .25  + (-!)(.  5) 


+ 1.8332  X (.75)  = .8333 
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step  (n+1 ) = 3 

i<2  = (-62/a2)  = -(-.0714/. 4286)  = .1666 


^30 

^ ^20 

= 1 

^31 

= 32^ 

+ k2 

322 

= -.8572 

+ .1666  X .1429  = 

-.8334 

^32 

= 322 

" ^2 

^1  = 

= .1429  + 

.1666  X (-.8572) 

= .0001 

®33 

= $2  ' 

*20  " 

.1666(1)  = .■ 

1666 

“3 

= a2  H 

H l<2 

62  = 

.4286  + , 

.1666  X (-.0714)  = 

.4167 

^3 

" ^30 

>"4^ 

^31 

"3  ^ "32 

"2  "33  >"1 

= 1 X 

(0)  • 

^ (-. 

,8334)(.25)  + (10‘^)(.5)  + 

(.1666)( 

°2 

= (93 

- ^2 

)/a3 

= (.25  - 

.8333)7.4167  = -1 

.3998 

^30  " '^'20  ^2  ^33  " -•''666  + (-1.3998)(.1666)  = -.3998 

^31  " ■^21  ^ ^2  ®32  " ^ (-1-3998)(10‘^)  = -1.0001 

^32  " ■''22  *^2  ^31  " ■'•8332  + (-1.3998) (-.8334)  = 2.9998 

^33  " ^2  ®30  " (-■'•3998)(1.)  = -1.3998 

^3  " ■''30  '"4  ^31  '"3  ^32  ’"2  ^33  '"l 

= -.3998  (0)  + (-1.001)(.25)  + (2.998)(.5)  + (-1 . 3998) ( . 75) 


= 2.000 
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step  (n+1) 

= 4 

: (-S3/a3)  = -(-.0834/. 41 67)  = 2.0001 

®40  " 

^30  " ^ 

^41 

^31  ^ *^3  ^33  ^ --8334  + (2. 0001 ) (.  1 666)  = -.8001 

®42  " 

’ ^32  4 ^32  " ■*■  (-20001  )(10'^)  = 10~^ 

®43  " 

®33  *^3  ®31  (-1 666)  + (. 20001 ) (-.8334)  = -.0002 

II 

<d 

IC3  83^  = (.2001)(1.)  = .2001 

^4  = 

“3  *^3  ^3  " (-.0834) (.2001  ) = .4000 

^4  = 

not  necessary 

^3  = 

(94  - Y3)/a4  = (0.  - (.2000))/. 4 = -.5000 

^40  " 

^30  ^3  ^44  " --3998  + (-.5)  x (.2001)  = -.4998 

^41  " 

^31  93  a^3  = -1.0001  + (-.5) (-.0002)  = -1. 

^42  " 

^32  * 93  842  = 2.9998  + (-.5)(.0001)  = 2.9997 

^43  " 

^33  93  84^  = -1.3998  + (-.5)(-.8001  ) = -.9998 

^44  " 

93  ^40  ~ (“-6)(1-)  - -.5000 

Solution: 

= -.4998 
= -1.0000 
f2  = 2.9997 
f3  = -.9998 
f4  = -.5000 
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The  result  can  be  checked  by  performing 
R . f = g 


1 

Lf) 

.5 

.25 

0 

.75 

1 

.75 

.5 

.25 

.5 

.75 

1 

.75 

.5 

.25 

.5 

.75 

1 

.75 

0 

.25 

.5 

.75 

1 

R 


-.4998 

0 

-1.000 

.25 

2.9998 

= 

.9999 

-.9999 

.2500 

-.5000 

0 

L 9 


The  starting  values  were 

9q  “ 0 , g-|  - . 25  , g2  ~ 1 » g^  ~ *25  , g^  = 


Therefore  f is  sought  after  solution. 


APPENDIX  B 

ORGANIGRAM  FOR  ITERATIVE  METHOD 
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A 


FF(L)  = DA' 

rA(l,L) 

> 

FF(L)  out 
FFF(L)  out  (F*^) 


* JM  is  an  odd  integer 


Figure  B.1  Organigram  for  the  iterative  technique. 


APPENDIX  C 

PROLATE  SPHEROIDAL  WAVE  FUNCTIONS 

The  prolate  spheroidal  wave  functions  are  scaled  versions  of  the 
angular  prolate  spheroidal  functions  which  are  solutions  of  the  differ- 
ential equation 

^ [(1  - t^)  ^ S(c,t)  + [X„„(c)  - ] S = 0 (C.l) 

where 

t in  [-1,+1] 

The  solution  of  the  first  kind  S||j^^(c,t)  is  finite,  analytic  in  [-1,+1] 
and  that  of  the  second  kind  Sj|j^^(c,t)  is  singular  at  t = ±1. 

C.l  Space-bandwidth  Parameter  c = 0 
The  differential  equation  becomes 

^ [(1  - t")  ^ S(o,t)]  + EXj^n^o)  - ] S(o,t)  = 0 (C.2) 

This  is  an  associated  Legendre  equation  whose  eigenvalues  are  given  by 

\<n  ' 

^ [(1  - t‘)  ^ S(o,t)]  + [n(n  + 1)  - = 0 (C.3) 

The  solution  is  a Legendre  polynomial  of  order  m and  degree  n 

Pj|(t)  = (1  - ^ P^(t)  (C.4) 

dt 
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P (t)  = — 


1 


2"  n! 


dt' 


(t^  - D" 


(C.5) 


C.2  Space-bandwidth  Parameter  c Not  Zero 
Equation  C.l  is  not  a Legendre  equation  and  has  a singularity  at 
infinity.  The  solution  in  this  case  can  be  found  in  the  form  [9]  of  an 
infinite  series: 


r=0,2,4,...  ^ 

IjSyS}. .* 


(C.6) 


The  summation  is  over  even  values  of  r for  (n  - m)  even  and  odd  values 


of  r for  (n  - m)  odd.  Where  the  expansion  coefficients  d|,^  are  related  by 
the  recurrence  relationship  given  by  Flammer  as  follows  (when  using  Equa- 
tion C.l  the  differential  equation  for  the  associated  Legendre  polynomials 
and  their  recurrence  relationship) 

( 2m  + r + 2 ) ( 2m  + r + 1 ) c ^ jfnn  , ^ , r ^ . n / . .-in 

(^m  t 2r  mk*  ir  * S)  * >•  + ' ) 

- Y (c)  + [2(m  + r)(m  + r + 1)  - 2m^  - 1]  c^  -|  ^mn 
(2m  + 2r  - l)(2m  + 2r  + 3) 


. r(r  - 1)  c^ .mn  ^ r, 

(2m  + 2r  - 3)(2m  + 2r  - 1 ) ^r-2^^^  ° ° 


(C.7) 


A normalization  condition  is  also  introduced 


(C.8) 


where 
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C<“) 


(n  ^ mil 


for  (n  - m)  even 


(C.9) 


and 


2^(0) 


n/n-m-^^,  /n+m+1 


for  (n  - m)  odd  (C.IO) 


2"  ( 


-)! 


)! 


C.3  Practical  Case,  m = 0 

A practical  case  of  interest  is  when  m = 0.  The  differential  equa- 
tion is  then 


^ [0  - t')  ^ Sg^Cc.t)]  + [Xq^(c)  - cH^]  S^^(c,t)  = 0 (c.n) 


t in  [-1,+1] 

The  solution  of  the  first  kind  is  S^^Cc.t)  corresponding  to  the  eigenvalue 
^on‘ 

C.3.1  Space-bandwidth  Parameter  c = 0 

It  [(’  - ^ + x„„(0)  S^„{o,t)  = 0 (C.12) 

The  solution  is  S^^Co.t)  = P^(t)  with  eigenvalues  given  by 

Xo„  ' + 1) 

and  the  normalization  condition  is 


14£ 

C.3.2  Space-bandwidth  Parameter  c Not  Zero 

The  solution  is  given  by  the  infinite  series  (Equation  C.6) 


S (c,t)  = I 
° r=0,2,4,, 

r=l ,3,5, , 


with  the  normalization  condition 


d°"(c)  P^(t) 


where 


and 


P („)  = 


for  n even 


The  SQ^(c,t)  function  has  exactly  n zeroes  in  the  interval 

t in  [-1,+1] 

Taking  advantage  of  this  fact,  the  eigenvalues  are  determined  in  a dif- 
ferent way  than  that  used  by  Flammer  in  the  form  of  an  infinite  series 


C.4  Eigenvalues  Normalized  to  One 

The  prolate  spheroidal  functions  are  also  the  solutions  of  the  homo- 
geneous Fredholm  integral  equation 


M{t)  = 


f+T 


(J)(x) 


-T 


sin  f2(t  - x) 
TT(t  - x) 


dx 


(C.13) 


with  the  following  properties: 
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C.4.1  Eigenvalues 

in  R , positive. 

1 > A-,  > X,  > ...  > X ».  0 

u I n 

as  n *■  oo 


C.4.2  Eigenfunctions 

To  a given  eigenvalue  X^  it  corresponds  an  eigenfunction  4>p(t)  such 

that 

*„(t)  *„(t)  dt  = (C.U) 

4 —00 

and 

= (-1)"  (|)^(t)  (C.15) 

And  they  are  also  orthogonal  in  another  interval  [8] 


(t)  <j)  (t)  dt  = X 6 
^n'  ' n mn 


(C.16) 


These  eigenfunctions  constitute  a complete  set.  Note:  the  space  bandwidth 

parameter  c is  defined  by 


c = 


T 


where  = Nyquist  limit. 


(C.17) 


APPENDIX  D 

DIGITAL  PROLATE  PROBLEMS 


In  this  section  a tentative  to  solve  the  prolate  spheroidal  wave 
equation  digitally  is  investigated.  As  mentioned  by  Frieden  Dl]  this 
method  has  never  been  approached  by  deconvolution  investigators. 


D.l  Digital  Deconvolution  and  Extrapolation 
As  noted  before,  the  angular  prolate  spheroidal  functions  are  the 
solutions  to  either  one  of  the  following  forms 


c'-  2t  (y  - c^t^) 

^n  ■ (1  - t2)  ^n  (1  - t2)  ^n 


S = 0 


or 


Vn(')  = 


+1 


_1  7r(t  - s)  ^n^^^ 


(D.l) 

(D.2) 


The  problem  is,  therefore,  in  a first  step  to  find  digitally  the  eigen- 
values from  the  integral  relationship  (Equation  D.2),  and  secondly  to 
digitally  find  the  eigenfunctions,  Sp(t)  for  t in  the  interval  [0,+t^]. 
Thirdly,  and  finally,  one  finds  the  prolate  spheroidal  functions  4>f,(t)  to 
deconvolve  and  extrapolate  a bandlimited  function  f(t). 

To  start,  let  us  consider  the  integral  relationship  (Equation  D.2) 
in  a matrix  formulation 

= " ^n  fn 

R = convolution  matrix 
To  alleviate  the  notation  one  represents 
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sin  c(t  - s)  ^ ^ 

nU  - s)'  “"y 

where 

t = (k  - l)h 
s = (j  - l)h 
and  h = 1/N 


For  t in  the  interval  [-1,+1] 
and  s in  the  interval  [-1,+1]  , 

The  matrix  R becomes: 


j 

s = -1 

0 -h  ...  -nh  -(n+l)h 


s = +1 

-(2n-l)h  -2nh 


k 


h 0 -h  -(n-l)h  -nh 


-(2n-2)h  -(2n-l)h 


nh  (n-l)h  (n-2)h  . . . +h  0 -h 


t=0 


(n+l)h  nh 


(n-1  )h 


2h  h 0 


(n+2)h  (n+l)h  nh  ...  2h  h 


-(n-l)h  -nh 
-h  ...  -(n-l)h 
0 -h  ...  -(n-2)h 


t=+l 


2nh  (2n-l)h  (2n-2)h  ... 


...  2h  h 0 


The  entries  of  the  matrix  are  (t  - s),  but  in  fact  they  are 
(sin  c(t  - s))/(7r(t  - s))  which  means  that  the  matrix  is  persymmetric 
and  Toeplitz. 


14<? 


D.2  Compact  Form  of  the  Convolution  Matrix 
Taking  advantage  of  the  symmetry  property  of  the  prolate  spheroidal 
functions  and  the  fact  that  the  convolution  matrix  is  persymmetric,  the 
latter  can  tie  reduced  to  half-size,  provided  that  the  number  of  mesh  points 
is  even,  and  still  be  symmetric  as  shown  below. 

Assuming  the  number  of  mesh  points  to  be  8,  the  integral  relationship 
(Equation  D.2)  can  be  written  as  follows: 


di  d2  dg  d^  dg  dg  d^  dg 

^4 

04 

d2  d.|  d2  dg  d^  dg  dg  dy 

'<=3 

^3 

dg  dg  d.|  dg  dg  d^  dg  dg 

02 

02 

^4  ^3  ^2  *^1  ^2  *^4  *^5 

01 

= X 

01 

^5  *^4  ^3  *^2  ^1  ^2  ^3  ^4 

±01 

n 

±01 

^6  ^5  ^4  ^3  ^2  ^1  ^2  ^3 

±02 

±02 

^7  ^6  ^5  ^4  ^3  ^2  ^1  ^2 

-‘*'3 

-'*'3 

^8  ^7  ^6  ^5  ^4  ^3  ^2  ^1 

±04 

±04 

The  matrix  ^(N,N)  is  collapsed  to  half-size  . 


■dr^8 

. d2^d? 

V''5' 

1 

-e- 

1 

•e- 

1 

d2±dy 

d3±d4 

“^3 

= X 

'*’3 

^3^^6 

d2±dg 

d,±d, 

*^2 

n 

<^2 

1 

Q. 

1+ 

Q. 

tn 

d3±d4 

d2±dg 

di±d2 

'f’l 

01 

The  new  matrix  is  symmetric.  The  + sign  is  for  even  order  of  the 

eigenvalue.  The  - sign  is  for  odd  order  of  the  eigenvalue. 
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A QL-algorithm  was  used  to  find  the  eigenvalues  of  this  collapsed 
matrix,  but  due  to  the  stiffness  of  the  problem  (i.e.,  the  eigenvalues 
were  widely  distributed)  it  was  impossible  to  obtain  the  eigenvalues  using 
this  digital  technique  and  we  believe  this  was  one  of  the  reasons  why 
other  researchers  did  not  approach  it  [20,21 ,22,23]. 

The  programs  involved  here  were: 


RECUR. FTN 

to  compute  the  vector  basis  for  the  convolution 

RECURS. FTN 

matrix 

to  find  the  inverse  matrix  R"^  using  a recursion 
algorithm  that  solves  a system  of  linear  equations 
R F = G, 

DIGPRO.FTN 

involving  EHOUSS.FTN  and  EQRT2S.FTN  for  the 
determination  of  the  eigenvalues  based  on  the 
QL-algorithm;  the  eigenvectors  were  obtained  with 
subroutine  EHOBKS.FTN. 

The  eigenvalues  may  be  obtained  from  = 1 down  to  A^  ~ 10"^  using  single 
precision  calculations  (four  byte  real  numbers),  and  down  to  A^  ~ 10'^®  in 
double  precision  (eight  byte  real  numbers).  Extended  precision  in  a main- 
frame computer  permitted  even  smaller  eigenvalues,  but  at  best  we  were  only 
able  to  obtain  the  first  twenty  eigenvalues  for  c = 50. 
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